Show code cell content
# Source the package setup script
source("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/scripts/00_setup_packages.R")
# Source the custom graphing functions
source("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/scripts/01_graphing_functions.R")
The C++ toolchain required for CmdStan is setup properly!
* Installing CmdStan from https://github.com/stan-dev/cmdstan/releases/download/v2.32.1/cmdstan-2.32.1.tar.gz
Warning message:
"An installation already exists at C:\Users\bmc82/.cmdstan/cmdstan-2.32.1. Please remove or rename the installation folder or set overwrite=TRUE."
Attaching package: 'dplyr'
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
Loading required package: usethis
Linking to libssh v0.11.0
Attaching package: 'kableExtra'
The following object is masked from 'package:dplyr':
group_rows
Warning message:
"package 'FSA' was built under R version 4.4.3"
Registered S3 methods overwritten by 'FSA':
method from
confint.boot car
hist.boot car
## FSA v0.9.6. See citation('FSA') if used in publication.
## Run fishR() for related website and fishR('IFAR') for related book.
Attaching package: 'plotly'
The following object is masked from 'package:ggplot2':
last_plot
The following object is masked from 'package:stats':
filter
The following object is masked from 'package:graphics':
layout
Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
Loading required package: rJava
Loading required package: leaps
Loading required package: sp
Attaching package: 'raster'
The following object is masked from 'package:plotly':
select
The following object is masked from 'package:dplyr':
select
Attaching package: 'recolorize'
The following object is masked from 'package:plotly':
add_image
Warning message:
"package 'webshot2' was built under R version 4.4.3"
This is bayesplot version 1.11.1
- Online documentation and vignettes at mc-stan.org/bayesplot
- bayesplot theme set to bayesplot::theme_default()
* Does _not_ affect other ggplot2 plots
* See ?bayesplot_theme_set for details on theme setting
Attaching package: 'bayesplot'
The following object is masked from 'package:plotrix':
plot_bg
This is posterior version 1.6.0
Attaching package: 'posterior'
The following object is masked from 'package:bayesplot':
rhat
The following objects are masked from 'package:raster':
%in%, match
The following objects are masked from 'package:stats':
mad, sd, var
The following objects are masked from 'package:base':
%in%, match
Attaching package: 'gridExtra'
The following object is masked from 'package:dplyr':
combine
Loading required package: Rcpp
Loading 'brms' package (version 2.22.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').
Attaching package: 'brms'
The following object is masked from 'package:bayesplot':
rhat
The following object is masked from 'package:stats':
ar
Attaching package: 'scales'
The following object is masked from 'package:plotrix':
rescale
── Attaching core tidyverse packages
✔ forcats 1.0.0 ✔ readr 2.1.5
✔ lubridate 1.9.4 ✔ stringr 1.5.1
✔ purrr 1.0.2
── Conflicts ───────────────────────
✖ readr::col_factor() masks scales::col_factor()
✖ gridExtra::combine() masks dplyr::combine()
✖ purrr::discard() masks scales::discard()
✖ raster::extract() masks tidyr::extract()
✖ plotly::filter() masks dplyr::filter(), stats::filter()
✖ kableExtra::group_rows() masks dplyr::group_rows()
✖ dplyr::lag() masks stats::lag()
✖ raster::select() masks plotly::select(), dplyr::select()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Attaching package: 'patchwork'
The following object is masked from 'package:raster':
area
Attaching package: 'cowplot'
The following object is masked from 'package:patchwork':
align_plots
The following object is masked from 'package:lubridate':
stamp
The following object is masked from 'package:ggpubr':
get_legend
This is loo version 2.8.0
- Online documentation and vignettes at mc-stan.org/loo
- As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session.
- Windows 10 users: loo may be very slow if 'mc.cores' is set in your .Rprofile file (see https://github.com/stan-dev/loo/issues/94).
Attaching package: 'extraDistr'
The following object is masked from 'package:purrr':
rdunif
The following objects are masked from 'package:brms':
ddirichlet, dfrechet, pfrechet, qfrechet, rdirichlet, rfrechet
Warning message:
"package 'leaflet' was built under R version 4.4.3"
Loading required package: StanHeaders
rstan version 2.32.6 (Stan version 2.32.2)
For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
change `threads_per_chain` option:
rstan_options(threads_per_chain = 1)
Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
Attaching package: 'rstan'
The following objects are masked from 'package:posterior':
ess_bulk, ess_tail
The following object is masked from 'package:raster':
extract
The following object is masked from 'package:tidyr':
extract
Warning message:
"package 'ggdist' was built under R version 4.4.3"
Attaching package: 'ggdist'
The following objects are masked from 'package:brms':
dstudent_t, pstudent_t, qstudent_t, rstudent_t
Attaching package: 'ape'
The following objects are masked from 'package:raster':
rotate, zoom
The following object is masked from 'package:glmulti':
consensus
The following object is masked from 'package:ggpubr':
rotate
The following object is masked from 'package:dplyr':
where
Bioconductor version '3.20' is out-of-date; the current release version '3.21'
is available with R version '4.5'; see https://bioconductor.org/install
Attaching package: 'BiocManager'
The following object is masked from 'package:devtools':
install
treeio v1.30.0 Learn more at https://yulab-smu.top/contribution-tree-data/
Please cite:
LG Wang, TTY Lam, S Xu, Z Dai, L Zhou, T Feng, P Guo, CW Dunn, BR
Jones, T Bradley, H Zhu, Y Guan, Y Jiang, G Yu. treeio: an R package
for phylogenetic tree input and output with richly annotated and
associated data. Molecular Biology and Evolution. 2020, 37(2):599-603.
doi: 10.1093/molbev/msz240
Attaching package: 'treeio'
The following object is masked from 'package:raster':
mask
Attaching package: 'TreeTools'
The following object is masked from 'package:treeio':
MRCA
ggtree v3.14.0 Learn more at https://yulab-smu.top/contribution-tree-data/
Please cite:
Guangchuang Yu. Data Integration, Manipulation and Visualization of
Phylogenetic Trees (1st edition). Chapman and Hall/CRC. 2022,
doi:10.1201/9781003279242, ISBN: 9781032233574
Attaching package: 'ggtree'
The following object is masked from 'package:TreeTools':
MRCA
The following object is masked from 'package:ape':
rotate
The following objects are masked from 'package:raster':
flip, rotate
The following object is masked from 'package:tidyr':
expand
The following object is masked from 'package:ggpubr':
rotate
Model 5#
Question#
Does camouflage strength (background matching) vary between males and females across microhabitats?
Objective#
Test for the effect of sex and microhabitat on background matching.
Method#
1. Load cleaned data.#
We start by loading the cleaned data from the “03_data_cleaning” pipeline. This data has already undergone transformations and contains relevant metrics for our models.
data_m5_clean <- read.csv("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/cleaned/data_m5_clean.csv")
2. Prepare data for modeling.#
Categorical predictors#
For categorical and binary predictors, R automatically dummy-codes the variables. For binary variables, such as Sex (0 = Female, 1 = Male), R sets the baseline to 0 unless specified otherwise. This ensures comparisons are made relative to the baseline group.
Continuous predictors#
We standardize continuous predictors, such as Fecundity, by centering and scaling them (dividing by two standard deviations). This improves interpretability, aligning their coefficients with those of binary predictors, as suggested by Gelman (2008).
# Convert categorical variables to factors
columns_to_convert_m5 <- c("Image", "Morph", "Sex", "Microhabitat_Association", "Viewpoint", "Microhabitat", "Background")
data_m5_clean <- data_m5_clean %>%
mutate(across(all_of(columns_to_convert_m5), as.factor))
# Set reference category for categorical predictors
data_m5_clean$Sex <- relevel(data_m5_clean$Sex, ref = "M")
data_m5_clean$Microhabitat <- relevel(data_m5_clean$Microhabitat, ref = "Hydroid")
data_m5_clean$Viewpoint <- relevel(data_m5_clean$Viewpoint, ref = "0")
# Convert continuous variables to numeric
columns_to_convert_m5 <- c("e_max_diff", "Filter_max_diff", "e_prop_diff", "R_c_diff", "G_c_diff", "B_c_diff")
data_m5_clean <- data_m5_clean %>%
mutate(across(all_of(columns_to_convert_m5), as.numeric))
3. Visualize response variable distributions#
To understand the distributions of our response variables, we plot density curves. This helps confirm whether a Gamma distribution is appropriate for modeling, as Gamma is suited for positively skewed, continuous data.
Show code cell content
p_response_m5_e_max_diff <- data_m5_clean %>%
ggplot(aes(x = e_max_diff)) +
geom_density(fill = "#69b3a2", color = "#e9ecef", alpha = 0.8) +
theme_bw() +
labs(title = "e_max_diff") +
theme(axis.title.x = element_blank())
p_response_m5_Filter_max_diff <- data_m5_clean %>%
ggplot(aes(x = Filter_max_diff)) +
geom_density(fill = "#69b3a2", color = "#e9ecef", alpha = 0.8) +
theme_bw() +
labs(title = "Filter_max_diff") +
theme(
axis.title.y = element_blank(),
axis.title.x = element_blank()
)
p_response_m5_e_prop_diff <- data_m5_clean %>%
ggplot(aes(x = e_prop_diff)) +
geom_density(fill = "#69b3a2", color = "#e9ecef", alpha = 0.8) +
theme_bw() +
labs(title = "e_prop_diff") +
theme(
axis.title.y = element_blank(),
axis.title.x = element_blank()
)
p_response_m5_R_c_diff <- data_m5_clean %>%
ggplot(aes(x = R_c_diff)) +
geom_density(fill = "#69b3a2", color = "#e9ecef", alpha = 0.8) +
theme_bw() +
labs(title = "R Reflectance diff") +
theme(axis.title.x = element_blank())
p_response_m5_G_c_diff <- data_m5_clean %>%
ggplot(aes(x = G_c_diff)) +
geom_density(fill = "#69b3a2", color = "#e9ecef", alpha = 0.8) +
theme_bw() +
labs(title = "G Reflectance diff") +
theme(
axis.title.y = element_blank(),
axis.title.x = element_blank()
)
p_response_m5_B_c_diff <- data_m5_clean %>%
ggplot(aes(x = B_c_diff)) +
geom_density(fill = "#69b3a2", color = "#e9ecef", alpha = 0.8) +
theme_bw() +
labs(title = "B Reflectance diff") +
theme(
axis.title.y = element_blank(),
axis.title.x = element_blank()
)
# Combine plots for easy visualization
plot_responses_m5 <- ggarrange(
p_response_m5_e_max_diff,
p_response_m5_Filter_max_diff,
p_response_m5_e_prop_diff,
p_response_m5_R_c_diff,
p_response_m5_G_c_diff,
p_response_m5_B_c_diff,
ncol = 3,
nrow = 2
)
annotate_figure(
plot_responses_m5,
top = text_grob("Ranges of response variables")
)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/plot_responses_m5.png", plot = plot_responses_m5, width = 6, height = 3, units = "in", dpi = 300)
Show code cell source
# Convert images to base64
plot_responses_m5 <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/plot_responses_m5.png")
# Create the HTML
html_background_matching <- paste0("
<style>
body, html {
margin: 0;
padding: 0;
/* If you want no horizontal scrollbar: */
overflow-x: hidden;
}
img {
max-width: 600px; /* ~6 inches at 100 dpi screen rendering */
width: 100%;
height: auto;
display: block;
margin-bottom: 20px;
border: 1px solid #ccc;
}
</style>
<img src='", plot_responses_m5, "' alt='Background Matching Metrics'>
")
# Display the HTML
IRdisplay::display_html(html_background_matching)
4. Justification for GLMs and Bayesian methods#
Why GLMs?#
Generalized Linear Models (GLMs) are used because our response variables are continuous, non-negative, and right-skewed, which violates the assumptions of linear regression. The Gamma distribution with a log link function allows us to model these variables appropriately.
Why Bayesian methods?#
Bayesian GLMs were chosen over frequentist approaches because:
Sparse data: Bayesian methods handle sparse datasets more robustly.
Priors: They allow us to incorporate prior knowledge, improving model performance.
Posterior distributions: Bayesian models provide posterior distributions, offering a full view of parameter uncertainty.
5. Define priors#
What are priors?#
In Bayesian statistics, priors represent our beliefs about parameter values before analyzing the data. These beliefs are mathematically expressed as probability distributions. Priors guide the model, especially when data are sparse or when the signal-to-noise ratio is low.
Why priors matter:#
Prevent overfitting: Priors discourage extreme parameter estimates unless strongly supported by the data.
Balance restrictiveness and flexibility: Weakly informative priors let the data dominate while providing reasonable bounds.
Leverage existing knowledge: Informative priors incorporate previous research or domain expertise, improving accuracy in well-studied systems.
How priors work in this analysis:#
We combine the priors (representing initial beliefs) with the likelihood of the observed data to compute posterior distributions, which reflect updated beliefs after observing the data.
The general formula is:
Model family and formula#
The response variables represent differences in color or pattern metrics between amphipod dorsal body traits and their microhabitats. These are continuous and positively skewed, making the Gamma distribution with a log link function suitable for modeling.
The Gamma distribution is parameterized as:
The mean \(\mu_i\) (on the log scale) is modeled as:
Where:
\(y_{i}\): Response variable for observation \(i\)
\(\mu_{i}\): Mean of the Gamma distribution for \(i\) (on the log scale)
\(\phi\): Shape parameter (controls variability)
\(\alpha\): Intercept, mean value of \(\log(\mu)\) when predictors are at baseline
\(beta\): Coefficient/slope for predictor \(x_{i}\), representing the effect of \(x_{i}\) on \(\log(\mu)\)
Chosen priors and rationale#
Intercept (\(\alpha\)): The intercept represents the mean of the response variable when all predictors are at their baseline values. Since our response variables are continuous and positively skewed, we use weakly informative priors for the intercept based on the log-transformed range of each variable. This ensures that the intercept (1) reflects realistic values for the response variable, and (2) allows for moderate variation around the mean without encouraging extreme values. For each response variable, we use the following priors:
e_max: \(N(0, 0.5)\)
\(N(0, 0.5)\) is a normal distribution with a mean of 0 and a standard deviation of \(\sqrt(0.5) = 0.707\). On the original scale, this translates to a mean of \(e^{0}=1\), and a standard deviation of \(e^{0.707}=2.03\). In a normal distribution, 95% of values fall within 2 standard deviations of the mean. Thus, most values will lie roughly between -1.414 and 1.414, on the log scale. On the original scale, the values lie roughly between \(e^{-1.414} \approx 0.243\) to \(e^{1.414} \approx 4.11\).
Filter_max: \(N(2.30, 0.5)\)
\(N(2.30, 0.5)\) is a normal distribution with a mean of 2.30 and a standard deviation of \(\sqrt(0.5) = 0.707\) On the original scale, this translates to a mean of \(e^{2.30}=10\), and a standard deviation of \(e^{0.707}=2.03\). In a normal distribution, 95% of values fall within 2 standard deviations of the mean. Thus, most values will lie roughly between 0.886 and 3.714, on the log scale. On the original scale, the values lie roughly between \(e^{0.886} \approx 2.43\) to \(e^{3.714} \approx 41\).
e_prop: \(N(-3.91, 0.5)\)
\(N(-3.91, 0.5)\) is a normal distribution with a mean of -3.91 and a standard deviation of \(\sqrt(0.5) = 0.707\). On the original scale, this translates to a mean of \(e^{-3.91}=0.02\), and a standard deviation of \(e^{0.707}=2.03\). In a normal distribution, 95% of values fall within 2 standard deviations of the mean. Thus, most values will lie roughly between -5.324 and -2.496, on the log scale. On the original scale, the values lie roughly between \(e^{-5.324} \approx 0.005\) to \(e^{-2.496} \approx 0.082\).
R, G, B: \(N(1.61, 0.5)\)
\(N(1.61, 0.5)\) is a normal distribution with a mean of 1.61 and a standard deviation of \(\sqrt(0.5) = 0.707\). On the original scale, this translates to a mean of \(e^{1.61}=5\), and a standard deviation of \(e^{0.707}=2.03\). In a normal distribution, 95% of values fall within 2 standard deviations of the mean. Thus, most values will lie roughly between 0.196 and 3.024, on the log scale. On the original scale, the values lie roughly between \(e^{0.196} \approx 1.22\) to \(e^{3.024} \approx 20.6\).
These priors ensure the intercept reflects realistic values for the response variables while allowing for moderate variability.
Slope (\(\beta\)): The slope reflects the rate of change in \(\log(\mu)\) for a one-unit change in the predictor. On the original scale, this translates to:
For example:
If \(\beta = 0.1\), a one-unit increase in the predictor scales \(\mu\) by \(e^{0.1} \approx 1.10\) (a 10% increase).
If \(\beta = -0.1\), a one-unit increase scales \(\mu\) by \(e^{-0.1} \approx 0.91\) (a 9% decrease).
We don’t know how our predictors will affect the response, so we use a weakly informative prior:
This prior assumes no effect of the predictor on average (\(e^{0} = 1\), and allows moderate positive or negative effects, spanning approximately \([-1, 1]\) on the log scale (\(e^{-1} \approx 0.37\) to \(e^{1} \approx 2.72\) on the original scale). So \(\mu\) of each response variable will not exceed a minimum of \(\mu * 0.37\) and maximum of \(\mu * 2.72\).
Random effects (\(\sigma\)): Random effects account for variability between groups that we are not explicitly testing. In this model, we do not have random effects.
Shape parameter (\(\phi\)): In a Gamma distribution, variability is tied to the shape parameter (\(\phi\)) that governs how the spread relates to the mean. Smaller \(\phi\) values indicate greater variability. We will set the shape parameter to a weakly informative \(\gamma(2, 0.1)\) which makes most probability mass concentrated near zero with a long tail so that the data will drive the distribution of variability values.
This ensures strictly positive values but keeps a wide range of plausible variability values.

Visualizing priors#
To validate these priors, we run models sampling only from the priors (sample_prior = “only”) and inspect their outputs to ensure they align with our expectations.
Note: We will run these models in RStudio to be consistent because the rstan package sometimes does not like Jupyter. The models were saved in RStudio and loaded below. We left the code chunks as comments for reference.
Show code cell content
# # Set seed for reproducibility
# set.seed(5678)
# # e_max_diff
# m5_priors_e_max_diff <- brm(
# bf(e_max_diff ~ 1 + Sex + Microhabitat + Viewpoint + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(0, 0.5), class = "Intercept"),
# prior(normal(0, 0.5), class = "b"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = "only",
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5_priors_e_max_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5_priors_e_max_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # Filter_max_diff
# m5_priors_Filter_max_diff <- brm(
# bf(Filter_max_diff ~ 1 + Sex + Microhabitat + Viewpoint + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(2.30, 0.5), class = "Intercept"),
# prior(normal(0, 0.5), class = "b"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = "only",
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5_priors_Filter_max_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5_priors_Filter_max_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # e_prop_diff
# m5_priors_e_prop_diff <- brm(
# bf(e_prop_diff ~ 1 + Sex + Microhabitat + Viewpoint + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(-3.91, 0.5), class = "Intercept"),
# prior(normal(0, 0.5), class = "b"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = "only",
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5_priors_e_prop_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5_priors_e_prop_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # R_c_diff
# m5_priors_R_c_diff <- brm(
# bf(R_c_diff ~ 1 + Sex + Microhabitat + Viewpoint + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(1.61, 0.5), class = "Intercept"),
# prior(normal(0, 0.5), class = "b"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = "only",
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5_priors_R_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5_priors_R_c_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # G_c_diff
# m5_priors_G_c_diff <- brm(
# bf(G_c_diff ~ 1 + Sex + Microhabitat + Viewpoint + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(1.61, 0.5), class = "Intercept"),
# prior(normal(0, 0.5), class = "b"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = "only",
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5_priors_G_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5_priors_G_c_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # B_c_diff
# m5_priors_B_c_diff <- brm(
# bf(B_c_diff ~ 1 + Sex + Microhabitat + Viewpoint + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(1.61, 0.5), class = "Intercept"),
# prior(normal(0, 0.5), class = "b"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = "only",
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5_priors_B_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5_priors_B_c_diff.rds")
Show code cell content
# Load models from R
m5_priors_e_max_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5_priors_e_max_diff.rds")
m5_priors_Filter_max_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5_priors_Filter_max_diff.rds")
m5_priors_e_prop_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5_priors_e_prop_diff.rds")
m5_priors_R_c_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5_priors_R_c_diff.rds")
m5_priors_G_c_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5_priors_G_c_diff.rds")
m5_priors_B_c_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5_priors_B_c_diff.rds")
# Extract results
prior_samples_m5_e_max_diff <- as_draws_df(m5_priors_e_max_diff)
prior_samples_m5_Filter_max_diff <- as_draws_df(m5_priors_Filter_max_diff)
prior_samples_m5_e_prop_diff <- as_draws_df(m5_priors_e_prop_diff)
prior_samples_m5_R_c_diff <- as_draws_df(m5_priors_R_c_diff)
prior_samples_m5_G_c_diff <- as_draws_df(m5_priors_G_c_diff)
prior_samples_m5_B_c_diff <- as_draws_df(m5_priors_B_c_diff)
Show code cell content
# List of datasets and labels
prior_samples <- list(
e_max_diff = prior_samples_m5_e_max_diff,
Filter_max_diff = prior_samples_m5_Filter_max_diff,
e_prop_diff = prior_samples_m5_e_prop_diff,
R_c_diff = prior_samples_m5_R_c_diff,
G_c_diff = prior_samples_m5_G_c_diff,
B_c_diff = prior_samples_m5_B_c_diff
)
# Custom labels for predictors
custom_labels_priors <- c(
"b_SexF" = "Females",
"b_MicrohabitatRed_Algae" = "Hydroids",
"b_MicrohabitatHydroid_Bryozoa" = "Hydroids with Bryozoans",
"b_Viewpoint10" = "Viewpoint 10",
"b_Viewpoint20" = "Viewpoint 20"
)
# Custom labels for intercept
custom_labels_priors_intercept <- c(
"b_Intercept" = "Intercept"
)
# Generate predictor plots
prior_plots <- lapply(
names(prior_samples),
function(label) {
generate_posterior_plot(
prior_samples[[label]],
regex_pars = c(
"b_SexF",
"b_MicrohabitatRed_Algae",
"b_MicrohabitatHydroid_Bryozoa",
"b_Viewpoint10",
"b_Viewpoint20"
),
x_range = c(-5, 5),
custom_labels = custom_labels_priors,
axis_title_y = label %in% c("e_max_diff", "R_c_diff")
)
}
)
names(prior_plots) <- names(prior_samples)
# Generate intercept plots
prior_plots_intercept <- lapply(
names(prior_samples),
function(label) {
generate_posterior_plot(
prior_samples[[label]],
regex_pars = c("b_Intercept"),
x_range = c(-5, 5),
custom_labels = custom_labels_priors_intercept,
axis_title_y = label %in% c("e_max_diff", "R_c_diff")
)
}
)
names(prior_plots_intercept) <- names(prior_samples)
# Add grey bars with labels for each plot
prior_plots_with_bars <- mapply(
function(plot, label) {
patchwork::wrap_elements(create_top_bar(label)) / plot +
patchwork::plot_layout(heights = c(0.2, 1))
},
prior_plots,
names(prior_plots),
SIMPLIFY = FALSE
)
prior_plots_intercept_with_bars <- mapply(
function(plot, label) {
patchwork::wrap_elements(create_top_bar(label)) / plot +
patchwork::plot_layout(heights = c(0.2, 1))
},
prior_plots_intercept,
names(prior_plots_intercept),
SIMPLIFY = FALSE
)
# Combine predictor plots into a 2x3 grid
plot_priors_m5 <- patchwork::wrap_plots(
prior_plots_with_bars[c("e_max_diff", "Filter_max_diff", "e_prop_diff", "R_c_diff", "G_c_diff", "B_c_diff")],
ncol = 3
) +
patchwork::plot_layout(guides = "collect")
# Combine intercept plots into a 2x3 grid
plot_priors_intercept_m5 <- patchwork::wrap_plots(
prior_plots_intercept_with_bars[c("e_max_diff", "Filter_max_diff", "e_prop_diff", "R_c_diff", "G_c_diff", "B_c_diff")],
ncol = 3
) +
patchwork::plot_layout(guides = "collect")
# Add a unified x-axis label for predictor plots
plot_priors_m5 <- plot_priors_m5 +
patchwork::plot_annotation(
caption = "Expected value of the response (log scale)",
theme = theme(plot.caption = element_text(hjust = 0.65, size = 10))
)
# Add a unified x-axis label for intercept plots
plot_priors_intercept_m5 <- plot_priors_intercept_m5 +
patchwork::plot_annotation(
caption = "Expected value of the response (log scale)",
theme = theme(plot.caption = element_text(hjust = 0.5, size = 10))
)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/plot_priors_m5.png", plot = plot_priors_m5, width = 6, height = 4.5, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/plot_priors_intercept_m5.png", plot = plot_priors_intercept_m5, width = 6, height = 4.5, units = "in", dpi = 300)
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Warning message:
"Using `size` aesthetic for lines
was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead."
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Warning message:
"Removed 1 row containing missing
values or values outside the scale
range (`geom_segment()`)."
Show code cell source
# Convert images to base64 (assuming these return base64 data URIs)
plot_priors_m5 <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/plot_priors_m5.png")
plot_priors_intercept_m5 <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/plot_priors_intercept_m5.png")
# Create the HTML
html_priors <- paste0("
<style>
.image-row {
display: flex;
gap: 20px;
justify-content: center;
align-items: flex-start;
}
.image-row img {
max-width: 100%;
height: auto;
border: 1px solid #ccc;
}
</style>
<div class='image-row'>
<img src='", plot_priors_m5, "' alt='Prior Plot'>
</div>
<div class='image-row'>
<img src='", plot_priors_intercept_m5, "' alt='Intercept Prior Plot'>
</div>
")
# Display the HTML
IRdisplay::display_html(html_priors)
6. Run final models#
Now that we have finalized the model parameters, we fit models using the actual data and compare them to null models to assess the significance of predictors.
Note: We will run these models in RStudio to be consistent because the rstan package sometimes does not like Jupyter. The models were saved in RStudio and loaded below. We left the code chunks as comments for reference.
# # Set seed for reproducibility
# set.seed(5678)
# # null model (for comparison)
# m5v0_e_max_diff <- brm(
# bf(e_max_diff ~ 1 + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(0, 0.5), class = "Intercept"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = TRUE,
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5v0_e_max_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v0_e_max_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # e_max_diff (Microhabitat)
# m5v1_e_max_diff <- brm(
# bf(e_max_diff ~ 1 + Sex + Microhabitat + Viewpoint + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(0, 0.5), class = "Intercept"),
# prior(normal(0, 0.5), class = "b"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = TRUE,
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5v1_e_max_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v1_e_max_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # e_max_diff (Microhabitat Association)
# m5v2_e_max_diff <- brm(
# bf(e_max_diff ~ 1 + Sex + Microhabitat_Association + Viewpoint + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(0, 0.5), class = "Intercept"),
# prior(normal(0, 0.5), class = "b"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = TRUE,
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5v2_e_max_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v2_e_max_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # null model (for comparison)
# m5v0_Filter_max_diff <- brm(
# bf(Filter_max_diff ~ 1 + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(2.30, 0.5), class = "Intercept"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = TRUE,
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical = FALSE),
# backend = "rstan"
# )
# saveRDS(m5v0_Filter_max_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v0_Filter_max_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # Filter_max_diff (Microhabitat)
# m5v1_Filter_max_diff <- brm(
# bf(Filter_max_diff ~ 1 + Sex + Microhabitat + Viewpoint + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(0, 0.5), class = "Intercept"),
# prior(normal(0, 0.5), class = "b"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = TRUE,
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5v1_Filter_max_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v1_Filter_max_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # Filter_max_diff (Microhabitat Association)
# m5v2_Filter_max_diff <- brm(
# bf(Filter_max_diff ~ 1 + Sex + Microhabitat_Association + Viewpoint + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(2.30, 0.5), class = "Intercept"),
# prior(normal(0, 0.5), class = "b"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = TRUE,
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical = FALSE),
# backend = "rstan"
# )
# saveRDS(m5v2_Filter_max_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v2_Filter_max_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # null model (for comparison)
# m5v0_e_prop_diff <- brm(
# bf(e_prop_diff ~ 1 + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(-3.91, 0.5), class = "Intercept"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = TRUE,
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5v0_e_prop_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v0_e_prop_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # e_prop_diff (Microhabitat)
# m5v1_e_prop_diff <- brm(
# bf(e_prop_diff ~ 1 + Sex + Microhabitat + Viewpoint + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(0, 0.5), class = "Intercept"),
# prior(normal(0, 0.5), class = "b"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = TRUE,
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5v1_e_prop_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v1_e_prop_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # e_prop_diff (Microhabitat Association)
# m5v2_e_prop_diff <- brm(
# bf(e_prop_diff ~ 1 + Sex + Microhabitat_Association + Viewpoint + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(-3.91, 0.5), class = "Intercept"),
# prior(normal(0, 0.5), class = "b"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = TRUE,
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5v2_e_prop_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v2_e_prop_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # null model (for comparison)
# m5v0_R_c_diff <- brm(
# bf(R_c_diff ~ 1 + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(1.61, 0.5), class = "Intercept"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = TRUE,
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5v0_R_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v0_R_c_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # R_c_diff (Microhabitat)
# m5v1_R_c_diff <- brm(
# bf(R_c_diff ~ 1 + Sex + Microhabitat + Viewpoint + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(0, 0.5), class = "Intercept"),
# prior(normal(0, 0.5), class = "b"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = TRUE,
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5v1_R_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v1_R_c_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # R_c_diff (Microhabitat Association)
# m5v2_R_c_diff <- brm(
# bf(R_c_diff ~ 1 + Sex + Microhabitat_Association + Viewpoint + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(1.61, 0.5), class = "Intercept"),
# prior(normal(0, 0.5), class = "b"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = TRUE,
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5v2_R_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v2_R_c_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # null model (for comparison)
# m5v0_G_c_diff <- brm(
# bf(G_c_diff ~ 1 + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(1.61, 0.5), class = "Intercept"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = TRUE,
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5v0_G_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v0_G_c_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # G_c_diff (Microhabitat)
# m5v1_G_c_diff <- brm(
# bf(G_c_diff ~ 1 + Sex + Microhabitat + Viewpoint + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(0, 0.5), class = "Intercept"),
# prior(normal(0, 0.5), class = "b"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = TRUE,
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5v1_G_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v1_G_c_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # G_c_diff (Microhabitat Association)
# m5v2_G_c_diff <- brm(
# bf(G_c_diff ~ 1 + Sex + Microhabitat_Association + Viewpoint + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(1.61, 0.5), class = "Intercept"),
# prior(normal(0, 0.5), class = "b"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = TRUE,
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5v2_G_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v2_G_c_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # null model (for comparison)
# m5v0_B_c_diff <- brm(
# bf(B_c_diff ~ 1 + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(1.61, 0.5), class = "Intercept"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = TRUE,
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5v0_B_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v0_B_c_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # B_c_diff (Microhabitat)
# m5v1_B_c_diff <- brm(
# bf(B_c_diff ~ 1 + Sex + Microhabitat + Viewpoint + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(0, 0.5), class = "Intercept"),
# prior(normal(0, 0.5), class = "b"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = TRUE,
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5v1_B_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v1_B_c_diff.rds")
# # Set seed for reproducibility
# set.seed(5678)
# # B_c_diff (Microhabitat Association)
# m5v2_B_c_diff <- brm(
# bf(B_c_diff ~ 1 + Sex + Microhabitat_Association + Viewpoint + (1|Image)),
# data = data_m5_clean,
# family = Gamma(link = "log"),
# prior = c(
# prior(normal(1.61, 0.5), class = "Intercept"),
# prior(normal(0, 0.5), class = "b"),
# prior(gamma(2, 0.1), class = "shape")
# ),
# sample_prior = TRUE,
# save_pars = save_pars(all = TRUE),
# control = list(adapt_delta = 0.99, max_treedepth = 12),
# iter = 15000,
# warmup = 5000,
# chains = 2,
# cores = parallel::detectCores(logical=FALSE),
# backend = "rstan"
# )
# saveRDS(m5v2_B_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v2_B_c_diff.rds")
Model comparison#
Make sure final models outperform the null model, and check to see which model performs best. Use LOO to compare models.
Note: We will added loo to these models in RStudio. The models were saved in RStudio and loaded below. We left the code chunks as comments for reference.
# # Add loo to models
# m5v0_e_max_diff <- add_criterion(m5v0_e_max_diff, "loo", moment_match = TRUE)
# m5v1_e_max_diff <- add_criterion(m5v1_e_max_diff, "loo", moment_match = TRUE)
# m5v2_e_max_diff <- add_criterion(m5v2_e_max_diff, "loo", moment_match = TRUE)
# m5v0_Filter_max_diff <- add_criterion(m5v0_Filter_max_diff, "loo", moment_match = TRUE)
# m5v1_Filter_max_diff <- add_criterion(m5v1_Filter_max_diff, "loo", moment_match = TRUE)
# m5v2_Filter_max_diff <- add_criterion(m5v2_Filter_max_diff, "loo", moment_match = TRUE)
# m5v0_e_prop_diff <- add_criterion(m5v0_e_prop_diff, "loo", moment_match = TRUE)
# m5v1_e_prop_diff <- add_criterion(m5v1_e_prop_diff, "loo", moment_match = TRUE)
# m5v2_e_prop_diff <- add_criterion(m5v2_e_prop_diff, "loo", moment_match = TRUE)
# m5v0_R_c_diff <- add_criterion(m5v0_R_c_diff, "loo", moment_match = TRUE)
# m5v1_R_c_diff <- add_criterion(m5v1_R_c_diff, "loo", moment_match = TRUE)
# m5v2_R_c_diff <- add_criterion(m5v2_R_c_diff, "loo", moment_match = TRUE)
# m5v0_G_c_diff <- add_criterion(m5v0_G_c_diff, "loo", moment_match = TRUE)
# m5v1_G_c_diff <- add_criterion(m5v1_G_c_diff, "loo", moment_match = TRUE)
# m5v2_G_c_diff <- add_criterion(m5v2_G_c_diff, "loo", moment_match = TRUE)
# m5v0_B_c_diff <- add_criterion(m5v0_B_c_diff, "loo", moment_match = TRUE)
# m5v1_B_c_diff <- add_criterion(m5v1_B_c_diff, "loo", moment_match = TRUE)
# m5v2_B_c_diff <- add_criterion(m5v2_B_c_diff, "loo", moment_match = TRUE)
# Save models with loo so you don't have to do this again
# saveRDS(m5v0_e_max_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v0_e_max_diff.rds")
# saveRDS(m5v1_e_max_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v1_e_max_diff.rds")
# saveRDS(m5v2_e_max_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v2_e_max_diff.rds")
# saveRDS(m5v0_Filter_max_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v0_Filter_max_diff.rds")
# saveRDS(m5v1_Filter_max_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v1_Filter_max_diff.rds")
# saveRDS(m5v2_Filter_max_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v2_Filter_max_diff.rds")
# saveRDS(m5v0_e_prop_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v0_e_prop_diff.rds")
# saveRDS(m5v1_e_prop_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v1_e_prop_diff.rds")
# saveRDS(m5v2_e_prop_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v2_e_prop_diff.rds")
# saveRDS(m5v0_R_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v0_R_c_diff.rds")
# saveRDS(m5v1_R_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v1_R_c_diff.rds")
# saveRDS(m5v2_R_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v2_R_c_diff.rds")
# saveRDS(m5v0_G_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v0_G_c_diff.rds")
# saveRDS(m5v1_G_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v1_G_c_diff.rds")
# saveRDS(m5v2_G_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v2_G_c_diff.rds")
# saveRDS(m5v0_B_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v0_B_c_diff.rds")
# saveRDS(m5v1_B_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v1_B_c_diff.rds")
# saveRDS(m5v2_B_c_diff, file = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v2_B_c_diff.rds")
Show code cell content
# Load models from R
m5v0_e_max_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v0_e_max_diff.rds")
m5v0_Filter_max_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v0_Filter_max_diff.rds")
m5v0_e_prop_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v0_e_prop_diff.rds")
m5v0_R_c_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v0_R_c_diff.rds")
m5v0_G_c_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v0_G_c_diff.rds")
m5v0_B_c_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v0_B_c_diff.rds")
m5v1_e_max_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v1_e_max_diff.rds")
m5v1_Filter_max_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v1_Filter_max_diff.rds")
m5v1_e_prop_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v1_e_prop_diff.rds")
m5v1_R_c_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v1_R_c_diff.rds")
m5v1_G_c_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v1_G_c_diff.rds")
m5v1_B_c_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v1_B_c_diff.rds")
m5v2_e_max_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v2_e_max_diff.rds")
m5v2_Filter_max_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v2_Filter_max_diff.rds")
m5v2_e_prop_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v2_e_prop_diff.rds")
m5v2_R_c_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v2_R_c_diff.rds")
m5v2_G_c_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v2_G_c_diff.rds")
m5v2_B_c_diff <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/models/m5v2_B_c_diff.rds")
# Compare all models
loo1_m5v1 <- loo_compare(m5v0_e_max_diff, m5v1_e_max_diff)
loo2_m5v1 <- loo_compare(m5v0_Filter_max_diff, m5v1_Filter_max_diff)
loo3_m5v1 <- loo_compare(m5v0_e_prop_diff, m5v1_e_prop_diff)
loo4_m5v1 <- loo_compare(m5v0_R_c_diff, m5v1_R_c_diff)
loo5_m5v1 <- loo_compare(m5v0_G_c_diff, m5v1_G_c_diff)
loo6_m5v1 <- loo_compare(m5v0_B_c_diff, m5v1_B_c_diff)
Show code cell content
# Convert to dataframe
df_loo1_m5v1 <- as.data.frame(loo1_m5v1) %>%
rownames_to_column(var = "Model")
df_loo2_m5v1 <- as.data.frame(loo2_m5v1) %>%
rownames_to_column(var = "Model")
df_loo3_m5v1 <- as.data.frame(loo3_m5v1) %>%
rownames_to_column(var = "Model")
df_loo4_m5v1 <- as.data.frame(loo4_m5v1) %>%
rownames_to_column(var = "Model")
df_loo5_m5v1 <- as.data.frame(loo5_m5v1) %>%
rownames_to_column(var = "Model")
df_loo6_m5v1 <- as.data.frame(loo6_m5v1) %>%
rownames_to_column(var = "Model")
# Save loo as tables
write.table(df_loo1_m5v1, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_loo_m5v1_e_max_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(df_loo2_m5v1, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_loo_m5v1_Filter_max_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(df_loo3_m5v1, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_loo_m5v1_e_prop_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(df_loo4_m5v1, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_loo_m5v1_R_c_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(df_loo5_m5v1, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_loo_m5v1_G_c_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(df_loo6_m5v1, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_loo_m5v1_B_c_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
Show code cell source
# Convert each data frame to a plain HTML table string
table_1_loo_m5v1 <- minimal_html_table(df_loo1_m5v1, caption = "LOO values m5v1 - Background Matching e_max")
table_2_loo_m5v1 <- minimal_html_table(df_loo2_m5v1, caption = "LOO values m5v1 - Background Matching Filter_max")
table_3_loo_m5v1 <- minimal_html_table(df_loo3_m5v1, caption = "LOO values m5v1 - Background Matching e_prop")
table_4_loo_m5v1 <- minimal_html_table(df_loo4_m5v1, caption = "LOO values m5v1 - Background Matching R Reflectance")
table_5_loo_m5v1 <- minimal_html_table(df_loo5_m5v1, caption = "LOO values m5v1 - Background Matching G Reflectance")
table_6_loo_m5v1 <- minimal_html_table(df_loo6_m5v1, caption = "LOO values m5v1 - Background Matching B Reflectance")
my_tabs_loo_m5v1 <- '
<style>
/* Basic container styling */
.tabs-container {
width: 100%;
margin: 1em 0;
}
/* Hide the radio inputs (we only show their labels as tabs) */
.tabs-container input[type="radio"] {
display: none;
}
/* The “tab-label” styling: looks like a tab */
.tab-label {
display: inline-block;
padding: 10px;
margin-right: 2px;
background: #eee;
border: 1px solid #ccc;
cursor: pointer;
border-bottom: none;
}
/* The active tab label */
.tab-label-active {
background: #fff;
}
/* The panel that holds table content */
.tab-content {
border: 1px solid #ccc;
padding: 10px;
display: none;
}
/* For each radio input, show its corresponding content when checked */
#tab1_loo_m5v1:checked ~ #content1_loo_m5v1,
#tab2_loo_m5v1:checked ~ #content2_loo_m5v1,
#tab3_loo_m5v1:checked ~ #content3_loo_m5v1,
#tab4_loo_m5v1:checked ~ #content4_loo_m5v1,
#tab5_loo_m5v1:checked ~ #content5_loo_m5v1,
#tab6_loo_m5v1:checked ~ #content6_loo_m5v1 {
display: block;
}
/* Also style the label of the checked radio as “active” using the :checked + label technique */
#tab1_loo_m5v1:checked + label[for="tab1_loo_m5v1"],
#tab2_loo_m5v1:checked + label[for="tab2_loo_m5v1"],
#tab3_loo_m5v1:checked + label[for="tab3_loo_m5v1"],
#tab4_loo_m5v1:checked + label[for="tab4_loo_m5v1"],
#tab5_loo_m5v1:checked + label[for="tab5_loo_m5v1"],
#tab6_loo_m5v1:checked + label[for="tab6_loo_m5v1"] {
background: #fff;
border-bottom: none;
}
</style>
<div class="tabs-container">
<!-- 1) Tab radio + label -->
<input type="radio" name="tabs_loo_m5v1" id="tab1_loo_m5v1" checked>
<label class="tab-label" for="tab1_loo_m5v1">Table 1</label>
<!-- 2) Tab radio + label -->
<input type="radio" name="tabs_loo_m5v1" id="tab2_loo_m5v1">
<label class="tab-label" for="tab2_loo_m5v1">Table 2</label>
<!-- 3) Tab radio + label -->
<input type="radio" name="tabs_loo_m5v1" id="tab3_loo_m5v1">
<label class="tab-label" for="tab3_loo_m5v1">Table 3</label>
<!-- 4) Tab radio + label -->
<input type="radio" name="tabs_loo_m5v1" id="tab4_loo_m5v1">
<label class="tab-label" for="tab4_loo_m5v1">Table 4</label>
<!-- 5) Tab radio + label -->
<input type="radio" name="tabs_loo_m5v1" id="tab5_loo_m5v1">
<label class="tab-label" for="tab5_loo_m5v1">Table 5</label>
<!-- 6) Tab radio + label -->
<input type="radio" name="tabs_loo_m5v1" id="tab6_loo_m5v1">
<label class="tab-label" for="tab6_loo_m5v1">Table 6</label>
<!-- Content for each tab -->
<div class="tab-content" id="content1_loo_m5v1">REPLACE_WITH_table_1_m5v1</div>
<div class="tab-content" id="content2_loo_m5v1">REPLACE_WITH_table_2_m5v1</div>
<div class="tab-content" id="content3_loo_m5v1">REPLACE_WITH_table_3_m5v1</div>
<div class="tab-content" id="content4_loo_m5v1">REPLACE_WITH_table_4_m5v1</div>
<div class="tab-content" id="content5_loo_m5v1">REPLACE_WITH_table_5_m5v1</div>
<div class="tab-content" id="content6_loo_m5v1">REPLACE_WITH_table_6_m5v1</div>
</div>
'
# Now do the replacements for each table
my_tabs_loo_m5v1 <- gsub("REPLACE_WITH_table_1_m5v1", table_1_loo_m5v1, my_tabs_loo_m5v1)
my_tabs_loo_m5v1 <- gsub("REPLACE_WITH_table_2_m5v1", table_2_loo_m5v1, my_tabs_loo_m5v1)
my_tabs_loo_m5v1 <- gsub("REPLACE_WITH_table_3_m5v1", table_3_loo_m5v1, my_tabs_loo_m5v1)
my_tabs_loo_m5v1 <- gsub("REPLACE_WITH_table_4_m5v1", table_4_loo_m5v1, my_tabs_loo_m5v1)
my_tabs_loo_m5v1 <- gsub("REPLACE_WITH_table_5_m5v1", table_5_loo_m5v1, my_tabs_loo_m5v1)
my_tabs_loo_m5v1 <- gsub("REPLACE_WITH_table_6_m5v1", table_6_loo_m5v1, my_tabs_loo_m5v1)
IRdisplay::display_html(my_tabs_loo_m5v1)
| Model | elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic |
|---|---|---|---|---|---|---|---|---|
| m5v1_e_max_diff | 0 | 0 | -2711.42003310893 | 133.243064576807 | 102.447609720779 | 1.6820073167562 | 5422.84006621785 | 266.486129153614 |
| m5v0_e_max_diff | -131.665157170932 | 14.651256827539 | -2843.08519027984 | 134.949856365264 | 99.1965204323952 | 1.70472416376797 | 5686.17038055969 | 269.899712730528 |
| Model | elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic |
|---|---|---|---|---|---|---|---|---|
| m5v1_Filter_max_diff | 0 | 0 | -46193.7932209322 | 130.789580277765 | 78.9269686550868 | 1.01233479315297 | 92387.5864418645 | 261.579160555531 |
| m5v0_Filter_max_diff | -8177.96045088377 | 100.569404835142 | -54371.7536718178 | 138.18103667034 | 98.6875547423149 | 1.76216847199398 | 108743.507343636 | 276.362073340681 |
| Model | elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic |
|---|---|---|---|---|---|---|---|---|
| m5v1_e_prop_diff | 0 | 0 | 60692.475827896 | 115.458210346294 | 91.2141461332322 | 1.25920210661846 | -121384.951655792 | 230.916420692589 |
| m5v0_e_prop_diff | -846.697287944139 | 39.1293967717956 | 59845.7785399521 | 115.459178890828 | 88.2051256335118 | 1.22733624096318 | -119691.557079904 | 230.918357781657 |
| Model | elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic |
|---|---|---|---|---|---|---|---|---|
| m5v1_R_c_diff | 0 | 0 | -46965.9229763443 | 136.876254185568 | 106.75664434762 | 2.18764806192843 | 93931.8459526887 | 273.752508371136 |
| m5v0_R_c_diff | -187.219569629834 | 20.1065004289648 | -47153.1425459742 | 139.733774235967 | 104.319466246147 | 2.36336282737204 | 94306.2850919484 | 279.467548471935 |
| Model | elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic |
|---|---|---|---|---|---|---|---|---|
| m5v1_G_c_diff | 0 | 0 | -46526.8133903886 | 135.374484552788 | 103.37674600417 | 1.95178736930457 | 93053.6267807771 | 270.748969105576 |
| m5v0_G_c_diff | -186.864783189826 | 18.5589362598883 | -46713.6781735782 | 137.496861912778 | 101.11115466969 | 2.05180225974255 | 93427.3563471564 | 274.993723825556 |
| Model | elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic |
|---|---|---|---|---|---|---|---|---|
| m5v1_B_c_diff | 0 | 0 | -42172.9130424138 | 150.496339733894 | 111.476507325309 | 2.02913676802372 | 84345.8260848276 | 300.992679467788 |
| m5v0_B_c_diff | -365.349125051042 | 26.3919836687448 | -42538.2621674647 | 152.86492648453 | 109.774306630455 | 2.28645015663594 | 85076.5243349293 | 305.72985296906 |
# Compare all models
loo1_m5v2 <- loo_compare(m5v0_e_max_diff, m5v2_e_max_diff)
loo2_m5v2 <- loo_compare(m5v0_Filter_max_diff, m5v2_Filter_max_diff)
loo3_m5v2 <- loo_compare(m5v0_e_prop_diff, m5v2_e_prop_diff)
loo4_m5v2 <- loo_compare(m5v0_R_c_diff, m5v2_R_c_diff)
loo5_m5v2 <- loo_compare(m5v0_G_c_diff, m5v2_G_c_diff)
loo6_m5v2 <- loo_compare(m5v0_B_c_diff, m5v2_B_c_diff)
Show code cell content
# Convert to dataframe
df_loo1_m5v2 <- as.data.frame(loo1_m5v2) %>%
rownames_to_column(var = "Model")
df_loo2_m5v2 <- as.data.frame(loo2_m5v2) %>%
rownames_to_column(var = "Model")
df_loo3_m5v2 <- as.data.frame(loo3_m5v2) %>%
rownames_to_column(var = "Model")
df_loo4_m5v2 <- as.data.frame(loo4_m5v2) %>%
rownames_to_column(var = "Model")
df_loo5_m5v2 <- as.data.frame(loo5_m5v2) %>%
rownames_to_column(var = "Model")
df_loo6_m5v2 <- as.data.frame(loo6_m5v2) %>%
rownames_to_column(var = "Model")
# Save loo as tables
write.table(df_loo1_m5v2, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_loo_m5v2_e_max_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(df_loo2_m5v2, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_loo_m5v2_Filter_max_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(df_loo3_m5v2, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_loo_m5v2_e_prop_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(df_loo4_m5v2, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_loo_m5v2_R_c_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(df_loo5_m5v2, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_loo_m5v2_G_c_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(df_loo6_m5v2, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_loo_m5v2_B_c_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
Show code cell source
# Convert each data frame to a plain HTML table string
table_1_loo_m5v2 <- minimal_html_table(df_loo1_m5v2, caption = "LOO values m5v2 - Background Matching e_max")
table_2_loo_m5v2 <- minimal_html_table(df_loo2_m5v2, caption = "LOO values m5v2 - Background Matching Filter_max")
table_3_loo_m5v2 <- minimal_html_table(df_loo3_m5v2, caption = "LOO values m5v2 - Background Matching e_prop")
table_4_loo_m5v2 <- minimal_html_table(df_loo4_m5v2, caption = "LOO values m5v2 - Background Matching R Reflectance")
table_5_loo_m5v2 <- minimal_html_table(df_loo5_m5v2, caption = "LOO values m5v2 - Background Matching G Reflectance")
table_6_loo_m5v2 <- minimal_html_table(df_loo6_m5v2, caption = "LOO values m5v2 - Background Matching B Reflectance")
my_tabs_loo_m5v2 <- '
<style>
/* Basic container styling */
.tabs-container {
width: 100%;
margin: 1em 0;
}
/* Hide the radio inputs (we only show their labels as tabs) */
.tabs-container input[type="radio"] {
display: none;
}
/* The “tab-label” styling: looks like a tab */
.tab-label {
display: inline-block;
padding: 10px;
margin-right: 2px;
background: #eee;
border: 1px solid #ccc;
cursor: pointer;
border-bottom: none;
}
/* The active tab label */
.tab-label-active {
background: #fff;
}
/* The panel that holds table content */
.tab-content {
border: 1px solid #ccc;
padding: 10px;
display: none;
}
/* For each radio input, show its corresponding content when checked */
#tab1_loo_m5v2:checked ~ #content1_loo_m5v2,
#tab2_loo_m5v2:checked ~ #content2_loo_m5v2,
#tab3_loo_m5v2:checked ~ #content3_loo_m5v2,
#tab4_loo_m5v2:checked ~ #content4_loo_m5v2,
#tab5_loo_m5v2:checked ~ #content5_loo_m5v2,
#tab6_loo_m5v2:checked ~ #content6_loo_m5v2 {
display: block;
}
/* Also style the label of the checked radio as “active” using the :checked + label technique */
#tab1_loo_m5v2:checked + label[for="tab1_loo_m5v2"],
#tab2_loo_m5v2:checked + label[for="tab2_loo_m5v2"],
#tab3_loo_m5v2:checked + label[for="tab3_loo_m5v2"],
#tab4_loo_m5v2:checked + label[for="tab4_loo_m5v2"],
#tab5_loo_m5v2:checked + label[for="tab5_loo_m5v2"],
#tab6_loo_m5v2:checked + label[for="tab6_loo_m5v2"] {
background: #fff;
border-bottom: none;
}
</style>
<div class="tabs-container">
<!-- 1) Tab radio + label -->
<input type="radio" name="tabs_loo_m5v2" id="tab1_loo_m5v2" checked>
<label class="tab-label" for="tab1_loo_m5v2">Table 1</label>
<!-- 2) Tab radio + label -->
<input type="radio" name="tabs_loo_m5v2" id="tab2_loo_m5v2">
<label class="tab-label" for="tab2_loo_m5v2">Table 2</label>
<!-- 3) Tab radio + label -->
<input type="radio" name="tabs_loo_m5v2" id="tab3_loo_m5v2">
<label class="tab-label" for="tab3_loo_m5v2">Table 3</label>
<!-- 4) Tab radio + label -->
<input type="radio" name="tabs_loo_m5v2" id="tab4_loo_m5v2">
<label class="tab-label" for="tab4_loo_m5v2">Table 4</label>
<!-- 5) Tab radio + label -->
<input type="radio" name="tabs_loo_m5v2" id="tab5_loo_m5v2">
<label class="tab-label" for="tab5_loo_m5v2">Table 5</label>
<!-- 6) Tab radio + label -->
<input type="radio" name="tabs_loo_m5v2" id="tab6_loo_m5v2">
<label class="tab-label" for="tab6_loo_m5v2">Table 6</label>
<!-- Content for each tab -->
<div class="tab-content" id="content1_loo_m5v2">REPLACE_WITH_table_1_m5v2</div>
<div class="tab-content" id="content2_loo_m5v2">REPLACE_WITH_table_2_m5v2</div>
<div class="tab-content" id="content3_loo_m5v2">REPLACE_WITH_table_3_m5v2</div>
<div class="tab-content" id="content4_loo_m5v2">REPLACE_WITH_table_4_m5v2</div>
<div class="tab-content" id="content5_loo_m5v2">REPLACE_WITH_table_5_m5v2</div>
<div class="tab-content" id="content6_loo_m5v2">REPLACE_WITH_table_6_m5v2</div>
</div>
'
# Now do the replacements for each table
my_tabs_loo_m5v2 <- gsub("REPLACE_WITH_table_1_m5v2", table_1_loo_m5v2, my_tabs_loo_m5v2)
my_tabs_loo_m5v2 <- gsub("REPLACE_WITH_table_2_m5v2", table_2_loo_m5v2, my_tabs_loo_m5v2)
my_tabs_loo_m5v2 <- gsub("REPLACE_WITH_table_3_m5v2", table_3_loo_m5v2, my_tabs_loo_m5v2)
my_tabs_loo_m5v2 <- gsub("REPLACE_WITH_table_4_m5v2", table_4_loo_m5v2, my_tabs_loo_m5v2)
my_tabs_loo_m5v2 <- gsub("REPLACE_WITH_table_5_m5v2", table_5_loo_m5v2, my_tabs_loo_m5v2)
my_tabs_loo_m5v2 <- gsub("REPLACE_WITH_table_6_m5v2", table_6_loo_m5v2, my_tabs_loo_m5v2)
IRdisplay::display_html(my_tabs_loo_m5v2)
| Model | elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic |
|---|---|---|---|---|---|---|---|---|
| m5v2_e_max_diff | 0 | 0 | -2784.90719838852 | 134.807674618251 | 103.765376297426 | 1.72928907243985 | 5569.81439677703 | 269.615349236502 |
| m5v0_e_max_diff | -58.1779918913145 | 12.4508493451589 | -2843.08519027984 | 134.949856365264 | 99.1965204323952 | 1.70472416376797 | 5686.17038055969 | 269.899712730528 |
| Model | elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic |
|---|---|---|---|---|---|---|---|---|
| m5v2_Filter_max_diff | 0 | 0 | -46230.31970728 | 131.339982035471 | 78.5087646907221 | 1.02585097702353 | 92460.63941456 | 262.679964070941 |
| m5v0_Filter_max_diff | -8141.43396453617 | 100.167963557035 | -54371.7536718178 | 138.18103667034 | 98.6875547423149 | 1.76216847199398 | 108743.507343636 | 276.362073340681 |
| Model | elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic |
|---|---|---|---|---|---|---|---|---|
| m5v2_e_prop_diff | 0 | 0 | 60505.6355211798 | 113.568589471179 | 89.6363722101689 | 1.20259889894136 | -121011.27104236 | 227.137178942357 |
| m5v0_e_prop_diff | -659.856981227993 | 36.8294708327522 | 59845.7785399521 | 115.459178890828 | 88.2051256335118 | 1.22733624096318 | -119691.557079904 | 230.918357781657 |
| Model | elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic |
|---|---|---|---|---|---|---|---|---|
| m5v2_R_c_diff | 0 | 0 | -47000.3745419683 | 137.976707932105 | 107.612676734088 | 2.23201230811384 | 94000.7490839366 | 275.953415864209 |
| m5v0_R_c_diff | -152.768004005989 | 19.4985479569406 | -47153.1425459742 | 139.733774235967 | 104.319466246147 | 2.36336282737204 | 94306.2850919484 | 279.467548471935 |
| Model | elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic |
|---|---|---|---|---|---|---|---|---|
| m5v2_G_c_diff | 0 | 0 | -46611.0492362456 | 136.385197014693 | 104.412959087539 | 1.98058530547232 | 93222.0984724913 | 272.770394029387 |
| m5v0_G_c_diff | -102.628937332746 | 15.9600085279904 | -46713.6781735782 | 137.496861912778 | 101.11115466969 | 2.05180225974255 | 93427.3563471564 | 274.993723825556 |
| Model | elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic |
|---|---|---|---|---|---|---|---|---|
| m5v2_B_c_diff | 0 | 0 | -42351.4943371316 | 151.266706364421 | 113.175143202646 | 2.07094954370071 | 84702.9886742632 | 302.533412728843 |
| m5v0_B_c_diff | -186.767830333441 | 22.8645602079982 | -42538.2621674647 | 152.86492648453 | 109.774306630455 | 2.28645015663594 | 85076.5243349293 | 305.72985296906 |
Visualize posteriors#
We extract the final model results and visualize the posteriors of our model parameters to get an idea of the significance of the results. This is what we did above when we were evaluating our priors.
Show code cell source
# Extract results
posterior_samples_m5v1_e_max_diff <- as_draws_df(m5v1_e_max_diff)
posterior_samples_m5v1_Filter_max_diff <- as_draws_df(m5v1_Filter_max_diff)
posterior_samples_m5v1_e_prop_diff <- as_draws_df(m5v1_e_prop_diff)
posterior_samples_m5v1_R_c_diff <- as_draws_df(m5v1_R_c_diff)
posterior_samples_m5v1_G_c_diff <- as_draws_df(m5v1_G_c_diff)
posterior_samples_m5v1_B_c_diff <- as_draws_df(m5v1_B_c_diff)
posterior_samples_m5v2_e_max_diff <- as_draws_df(m5v2_e_max_diff)
posterior_samples_m5v2_Filter_max_diff <- as_draws_df(m5v2_Filter_max_diff)
posterior_samples_m5v2_e_prop_diff <- as_draws_df(m5v2_e_prop_diff)
posterior_samples_m5v2_R_c_diff <- as_draws_df(m5v2_R_c_diff)
posterior_samples_m5v2_G_c_diff <- as_draws_df(m5v2_G_c_diff)
posterior_samples_m5v2_B_c_diff <- as_draws_df(m5v2_B_c_diff)
Show code cell content
# List of datasets and labels
posterior_samples <- list(
e_max_diff = posterior_samples_m5v1_e_max_diff,
Filter_max_diff = posterior_samples_m5v1_Filter_max_diff,
e_prop_diff = posterior_samples_m5v1_e_prop_diff,
R_c_diff = posterior_samples_m5v1_R_c_diff,
G_c_diff = posterior_samples_m5v1_G_c_diff,
B_c_diff = posterior_samples_m5v1_B_c_diff
)
# Baseline category data
baseline_data <- tibble(
parameter = c("Males", "Hydroids", "No correction"),
mean = 0, # Centered at 0
ci_low = -0.2, # Example CI range
ci_high = 0.2
)
# Order categories
parameter_order <- c(
"Males",
"b_SexF",
"Hydroids",
"b_MicrohabitatRed_Algae",
"b_MicrohabitatHydroid_Bryozoa",
"No correction",
"b_Viewpoint10",
"b_Viewpoint20"
)
#Custom labels
custom_labels_posteriors <- c(
"Males" = "Males",
"b_SexF" = "Females",
"Hydroids" = "Hydroids",
"b_MicrohabitatRed_Algae" = "Red Algae",
"b_MicrohabitatHydroid_Bryozoa" = "Hydroids with Bryozoa",
"No correction" = "No correction",
"b_Viewpoint10" = "Viewpoint 10",
"b_Viewpoint20" = "Viewpoint 20"
)
# Convert the y-axis parameters to factors for consistent alignment
baseline_data$parameter <- factor(baseline_data$parameter, levels = parameter_order)
# Generate plots for predictors
posterior_plots <- lapply(
names(posterior_samples),
function(label) {
generate_posterior_plot(
posterior_samples[[label]],
regex_pars = c(
"b_SexF",
"b_MicrohabitatRed_Algae",
"b_MicrohabitatHydroid_Bryozoa",
"b_Viewpoint10",
"b_Viewpoint20"
),
x_range = c(-0.35, 0.35),
custom_labels = custom_labels_posteriors,
axis_title_y = label %in% c("e_max_diff", "R_c_diff"),
axis_text_x = label %in% c("R_c_diff", "G_c_diff", "B_c_diff") # Show x-axis text only for second-row labels
) +
geom_point(
data = baseline_data,
aes(x = mean, y = parameter),
inherit.aes = FALSE,
color = "dodgerblue4",
size = 2
)
}
)
names(posterior_plots) <- names(posterior_samples)
# Add grey bars with labels for each plot
posterior_plots_with_bars <- mapply(
function(plot, label) {
patchwork::wrap_elements(create_top_bar(label)) / plot +
patchwork::plot_layout(heights = c(0.2, 1))
},
posterior_plots,
names(posterior_plots),
SIMPLIFY = FALSE
)
# Combine predictor plots into a 2x3 grid
plot_posteriors_m5v1 <- patchwork::wrap_plots(
posterior_plots_with_bars[c("e_max_diff", "Filter_max_diff", "e_prop_diff", "R_c_diff", "G_c_diff", "B_c_diff")],
ncol = 3
) +
patchwork::plot_layout(guides = "collect")
# Add a unified x-axis label for predictor plots
plot_posteriors_m5v1 <- plot_posteriors_m5v1 +
patchwork::plot_annotation(
caption = "Expected value of the response (log scale)",
theme = theme(plot.caption = element_text(hjust = 0.55, size = 10))
)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/plot_posteriors_m5v1.png", plot = plot_posteriors_m5v1, width = 6, height = 4.5, units = "in", dpi = 300)
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Warning message:
"Removed 2 rows containing missing
values or values outside the scale
range (`geom_segment()`)."
Warning message:
"Removed 1 row containing missing
values or values outside the scale
range (`geom_segment()`)."
Show code cell source
# Convert images to base64
plot_posteriors_m5v1 <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/plot_posteriors_m5v1.png")
# Create the HTML
html_posteriors_m5v1 <- paste0("
<style>
.image-row {
display: flex;
gap: 20px;
justify-content: center;
align-items: flex-start;
}
.image-row img {
max-width: 100%;
height: auto;
border: 1px solid #ccc;
}
</style>
<img src='", plot_posteriors_m5v1, "' alt='Posterior Plot (Model 2)'>
")
# Display the HTML
IRdisplay::display_html(html_posteriors_m5v1)
Show code cell content
# List of datasets and labels
posterior_samples <- list(
e_max_diff = posterior_samples_m5v2_e_max_diff,
Filter_max_diff = posterior_samples_m5v2_Filter_max_diff,
e_prop_diff = posterior_samples_m5v2_e_prop_diff,
R_c_diff = posterior_samples_m5v2_R_c_diff,
G_c_diff = posterior_samples_m5v2_G_c_diff,
B_c_diff = posterior_samples_m5v2_B_c_diff
)
# Baseline category data
baseline_data <- tibble(
parameter = c("Males", "Microhabitat_Association Not Present", "No correction"),
mean = 0, # Centered at 0
ci_low = -0.2, # Example CI range
ci_high = 0.2
)
# Order categories
parameter_order <- c(
"Males",
"b_SexF",
"Hydroids",
"Microhabitat_Association Not Present",
"b_Microhabitat_AssociationPresent",
"No correction",
"b_Viewpoint10",
"b_Viewpoint20"
)
#Custom labels
custom_labels_posteriors <- c(
"Males" = "Males",
"b_SexF" = "Females",
"Hydroids" = "Hydroids",
"Microhabitat_Association Not Present" = "Unassociated",
"b_Microhabitat_AssociationPresent" = "Associated",
"No correction" = "No correction",
"b_Viewpoint10" = "Viewpoint 10",
"b_Viewpoint20" = "Viewpoint 20"
)
# Convert the y-axis parameters to factors for consistent alignment
baseline_data$parameter <- factor(baseline_data$parameter, levels = parameter_order)
# Generate plots for predictors
posterior_plots <- lapply(
names(posterior_samples),
function(label) {
generate_posterior_plot(
posterior_samples[[label]],
regex_pars = c(
"b_SexF",
"b_Microhabitat_AssociationPresent",
"b_Viewpoint10",
"b_Viewpoint20"
),
x_range = c(-0.35, 0.35),
custom_labels = custom_labels_posteriors,
axis_title_y = label %in% c("e_max_diff", "R_c_diff"),
axis_text_x = label %in% c("R_c_diff", "G_c_diff", "B_c_diff") # Show x-axis text only for second-row labels
) +
geom_point(
data = baseline_data,
aes(x = mean, y = parameter),
inherit.aes = FALSE,
color = "dodgerblue4",
size = 2
)
}
)
names(posterior_plots) <- names(posterior_samples)
# Add grey bars with labels for each plot
posterior_plots_with_bars <- mapply(
function(plot, label) {
patchwork::wrap_elements(create_top_bar(label)) / plot +
patchwork::plot_layout(heights = c(0.2, 1))
},
posterior_plots,
names(posterior_plots),
SIMPLIFY = FALSE
)
# Combine predictor plots into a 2x3 grid
plot_posteriors_m5v2 <- patchwork::wrap_plots(
posterior_plots_with_bars[c("e_max_diff", "Filter_max_diff", "e_prop_diff", "R_c_diff", "G_c_diff", "B_c_diff")],
ncol = 3
) +
patchwork::plot_layout(guides = "collect")
# Add a unified x-axis label for predictor plots
plot_posteriors_m5v2 <- plot_posteriors_m5v2 +
patchwork::plot_annotation(
caption = "Expected value of the response (log scale)",
theme = theme(plot.caption = element_text(hjust = 0.55, size = 10))
)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/plot_posteriors_m5v2.png", plot = plot_posteriors_m5v2, width = 6, height = 4.5, units = "in", dpi = 300)
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which
will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which
will replace the existing scale.
Warning message:
"Removed 2 rows containing missing
values or values outside the scale
range (`geom_segment()`)."
Show code cell source
# Convert images to base64
plot_posteriors_m5v2 <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/plot_posteriors_m5v2.png")
# Create the HTML
html_posteriors_m5v2 <- paste0("
<style>
.image-row {
display: flex;
gap: 20px;
justify-content: center;
align-items: flex-start;
}
.image-row img {
max-width: 100%;
height: auto;
border: 1px solid #ccc;
}
</style>
<img src='", plot_posteriors_m5v2, "' alt='Posterior Plot (Model 2)'>
")
# Display the HTML
IRdisplay::display_html(html_posteriors_m5v2)
The graph of the posteriors gives us an idea of the significance of each predictor. We need to follow up with an evaluation of the model performance before we can trust these results.
7. Evaluate model performance#
Trace plots#
Visualize parameter sampling across iterations to confirm convergence. Each chain should wander around the same mean value without any strong upward or downward trends. “fuzzy caterpillar” or “horizontal band.”
Show code cell content
# Generate trace plots for all models
trace_plots_m5v1_e_max_diff <- generate_trace_plot(
model = m5v1_e_max_diff,
regex_pars = c("b_SexF", "b_MicrohabitatRed_Algae", "b_MicrohabitatHydroid_Bryozoa", "b_Viewpoint10", "b_Viewpoint20"),
plot_title = "e_max Model",
axis_title_y = TRUE,
axis_text_x = TRUE,
axis_title_size = 10,
axis_text_size = 8
)
trace_plots_m5v1_Filter_max_diff <- generate_trace_plot(
model = m5v1_Filter_max_diff,
regex_pars = c("b_SexF", "b_MicrohabitatRed_Algae", "b_MicrohabitatHydroid_Bryozoa", "b_Viewpoint10", "b_Viewpoint20"),
plot_title = "Filter_max Model",
axis_title_y = TRUE,
axis_text_x = TRUE,
axis_title_size = 10,
axis_text_size = 8
)
trace_plots_m5v1_e_prop_diff <- generate_trace_plot(
model = m5v1_e_prop_diff,
regex_pars = c("b_SexF", "b_MicrohabitatRed_Algae", "b_MicrohabitatHydroid_Bryozoa", "b_Viewpoint10", "b_Viewpoint20"),
plot_title = "e_prop Model",
axis_title_y = TRUE,
axis_text_x = TRUE,
axis_title_size = 10,
axis_text_size = 8
)
trace_plots_m5v1_R_c_diff <- generate_trace_plot(
model = m5v1_R_c_diff,
regex_pars = c("b_SexF", "b_MicrohabitatRed_Algae", "b_MicrohabitatHydroid_Bryozoa", "b_Viewpoint10", "b_Viewpoint20"),
plot_title = "R Reflectance Model",
axis_title_y = TRUE,
axis_text_x = TRUE,
axis_title_size = 10,
axis_text_size = 8
)
trace_plots_m5v1_G_c_diff <- generate_trace_plot(
model = m5v1_G_c_diff,
regex_pars = c("b_SexF", "b_MicrohabitatRed_Algae", "b_MicrohabitatHydroid_Bryozoa", "b_Viewpoint10", "b_Viewpoint20"),
plot_title = "G Reflectance Model",
axis_title_y = TRUE,
axis_text_x = TRUE,
axis_title_size = 10,
axis_text_size = 8
)
trace_plots_m5v1_B_c_diff <- generate_trace_plot(
model = m5v1_B_c_diff,
regex_pars = c("b_SexF", "b_MicrohabitatRed_Algae", "b_MicrohabitatHydroid_Bryozoa", "b_Viewpoint10", "b_Viewpoint20"),
plot_title = "B Reflectance Model",
axis_title_y = TRUE,
axis_text_x = TRUE,
axis_title_size = 10,
axis_text_size = 8
)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v1_e_max_diff.png", plot = trace_plots_m5v1_e_max_diff, width = 6, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v1_Filter_max_diff.png", plot = trace_plots_m5v1_Filter_max_diff, width = 6, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v1_e_prop_diff.png", plot = trace_plots_m5v1_e_prop_diff, width = 6, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v1_R_c_diff.png", plot = trace_plots_m5v1_R_c_diff, width = 6, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v1_G_c_diff.png", plot = trace_plots_m5v1_G_c_diff, width = 6, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v1_B_c_diff.png", plot = trace_plots_m5v1_B_c_diff, width = 6, height = 3, units = "in", dpi = 300)
Show code cell source
# Convert images to base64 (if not already done)
trace_plots_m5v1_e_max_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v1_e_max_diff.png")
trace_plots_m5v1_Filter_max_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v1_Filter_max_diff.png")
trace_plots_m5v1_e_prop_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v1_e_prop_diff.png")
trace_plots_m5v1_R_c_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v1_R_c_diff.png")
trace_plots_m5v1_G_c_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v1_G_c_diff.png")
trace_plots_m5v1_B_c_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v1_B_c_diff.png")
# Create the HTML
html_trace_plots_m5v1 <- paste0("
<style>
.grid-container {
display: grid;
grid-template-columns: repeat(2, 1fr); /* 2 columns per row */
gap: 3px;
padding: 3px;
justify-items: center;
}
.grid-container img {
max-width: 600px;
width: 100%;
height: auto;
border: 1px solid #ccc;
}
</style>
<div class='grid-container'>
<img src='", trace_plots_m5v1_e_max_diff, "' alt='e_max Plot (Model 1)'>
<img src='", trace_plots_m5v1_Filter_max_diff, "' alt='Filter_max Plot (Model 1)'>
<img src='", trace_plots_m5v1_e_prop_diff, "' alt='e_prop Plot (Model 1)'>
<img src='", trace_plots_m5v1_R_c_diff, "' alt='R Reflectance Plot (Model 1)'>
<img src='", trace_plots_m5v1_G_c_diff, "' alt='G Reflectance Plot (Model 1)'>
<img src='", trace_plots_m5v1_B_c_diff, "' alt='B Reflectance Plot (Model 1)'>
</div>
")
# Display the HTML
IRdisplay::display_html(html_trace_plots_m5v1)
Show code cell content
# Generate trace plots for all models
trace_plots_m5v2_e_max_diff <- generate_trace_plot(
model = m5v2_e_max_diff,
regex_pars = c("b_SexF", "b_Microhabitat_AssociationPresent", "b_Viewpoint10", "b_Viewpoint20"),
plot_title = "e_max Model",
axis_title_y = TRUE,
axis_text_x = TRUE,
axis_title_size = 10,
axis_text_size = 8
)
trace_plots_m5v2_Filter_max_diff <- generate_trace_plot(
model = m5v2_Filter_max_diff,
regex_pars = c("b_SexF", "b_Microhabitat_AssociationPresent", "b_Viewpoint10", "b_Viewpoint20"),
plot_title = "Filter_max Model",
axis_title_y = TRUE,
axis_text_x = TRUE,
axis_title_size = 10,
axis_text_size = 8
)
trace_plots_m5v2_e_prop_diff <- generate_trace_plot(
model = m5v2_e_prop_diff,
regex_pars = c("b_SexF", "b_Microhabitat_AssociationPresent", "b_Viewpoint10", "b_Viewpoint20"),
plot_title = "e_prop Model",
axis_title_y = TRUE,
axis_text_x = TRUE,
axis_title_size = 10,
axis_text_size = 8
)
trace_plots_m5v2_R_c_diff <- generate_trace_plot(
model = m5v2_R_c_diff,
regex_pars = c("b_SexF", "b_Microhabitat_AssociationPresent", "b_Viewpoint10", "b_Viewpoint20"),
plot_title = "R Reflectance Model",
axis_title_y = TRUE,
axis_text_x = TRUE,
axis_title_size = 10,
axis_text_size = 8
)
trace_plots_m5v2_G_c_diff <- generate_trace_plot(
model = m5v2_G_c_diff,
regex_pars = c("b_SexF", "b_Microhabitat_AssociationPresent", "b_Viewpoint10", "b_Viewpoint20"),
plot_title = "G Reflectance Model",
axis_title_y = TRUE,
axis_text_x = TRUE,
axis_title_size = 10,
axis_text_size = 8
)
trace_plots_m5v2_B_c_diff <- generate_trace_plot(
model = m5v2_B_c_diff,
regex_pars = c("b_SexF", "b_Microhabitat_AssociationPresent", "b_Viewpoint10", "b_Viewpoint20"),
plot_title = "B Reflectance Model",
axis_title_y = TRUE,
axis_text_x = TRUE,
axis_title_size = 10,
axis_text_size = 8
)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v2_e_max_diff.png", plot = trace_plots_m5v2_e_max_diff, width = 6, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v2_Filter_max_diff.png", plot = trace_plots_m5v2_Filter_max_diff, width = 6, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v2_e_prop_diff.png", plot = trace_plots_m5v2_e_prop_diff, width = 6, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v2_R_c_diff.png", plot = trace_plots_m5v2_R_c_diff, width = 6, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v2_G_c_diff.png", plot = trace_plots_m5v2_G_c_diff, width = 6, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v2_B_c_diff.png", plot = trace_plots_m5v2_B_c_diff, width = 6, height = 3, units = "in", dpi = 300)
Show code cell source
# Convert images to base64 (if not already done)
trace_plots_m5v2_e_max_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v2_e_max_diff.png")
trace_plots_m5v2_Filter_max_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v2_Filter_max_diff.png")
trace_plots_m5v2_e_prop_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v2_e_prop_diff.png")
trace_plots_m5v2_R_c_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v2_R_c_diff.png")
trace_plots_m5v2_G_c_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v2_G_c_diff.png")
trace_plots_m5v2_B_c_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/trace_plots_m5v2_B_c_diff.png")
# Create the HTML
html_trace_plots_m5v2 <- paste0("
<style>
.grid-container {
display: grid;
grid-template-columns: repeat(2, 1fr); /* 2 columns per row */
gap: 3px;
padding: 3px;
justify-items: center;
}
.grid-container img {
max-width: 600px;
width: 100%;
height: auto;
border: 1px solid #ccc;
}
</style>
<div class='grid-container'>
<img src='", trace_plots_m5v2_e_max_diff, "' alt='e_max Plot (Model 2)'>
<img src='", trace_plots_m5v2_Filter_max_diff, "' alt='Filter_max Plot (Model 2)'>
<img src='", trace_plots_m5v2_e_prop_diff, "' alt='e_prop Plot (Model 2)'>
<img src='", trace_plots_m5v2_R_c_diff, "' alt='R Reflectance Plot (Model 2)'>
<img src='", trace_plots_m5v2_G_c_diff, "' alt='G Reflectance Plot (Model 2)'>
<img src='", trace_plots_m5v2_B_c_diff, "' alt='B Reflectance Plot (Model 2)'>
</div>
")
# Display the HTML
IRdisplay::display_html(html_trace_plots_m5v2)
Posterior predictive checks#
Simulate data based on the models and compare to observed data to verify goodness-of-fit. We do this using “pp_check”. The observed data (black line/dots) should sit comfortably within the distribution of simulated data (colored areas or lines)
Show code cell content
# Generate posterior predictive check plots for all models
pp_check_plots_m5v1_e_max_diff <- generate_pp_check(
model = m5v1_e_max_diff,
nreps = 100,
axis_title_y = TRUE,
y_label = "Density",
plot_title = "e_max Background Matching Model",
axis_title_size = 10, # custom title size
axis_text_size = 8 # custom tick label size
)
pp_check_plots_m5v1_Filter_max_diff <- generate_pp_check(
model = m5v1_Filter_max_diff,
nreps = 100,
axis_title_y = TRUE,
y_label = "Density",
plot_title = "Filter_max Background Matching Model",
axis_title_size = 10, # custom title size
axis_text_size = 8 # custom tick label size
)
# Generate posterior predictive check plots for all models
pp_check_plots_m5v1_e_prop_diff <- generate_pp_check(
model = m5v1_e_prop_diff,
nreps = 100,
axis_title_y = TRUE,
y_label = "Density",
plot_title = "e_prop Background Matching Model",
axis_title_size = 10, # custom title size
axis_text_size = 8 # custom tick label size
)
pp_check_plots_m5v1_R_c_diff <- generate_pp_check(
model = m5v1_R_c_diff,
nreps = 100,
axis_title_y = TRUE,
y_label = "Density",
plot_title = "R Reflectance Background Matching Model",
axis_title_size = 10, # custom title size
axis_text_size = 8 # custom tick label size
)
# Generate posterior predictive check plots for all models
pp_check_plots_m5v1_G_c_diff <- generate_pp_check(
model = m5v1_G_c_diff,
nreps = 100,
axis_title_y = TRUE,
y_label = "Density",
plot_title = "G Reflectance Background Matching Model",
axis_title_size = 10, # custom title size
axis_text_size = 8 # custom tick label size
)
pp_check_plots_m5v1_B_c_diff <- generate_pp_check(
model = m5v1_B_c_diff,
nreps = 100,
axis_title_y = TRUE,
y_label = "Density",
plot_title = "B Reflectance Background Matching Model",
axis_title_size = 10, # custom title size
axis_text_size = 8 # custom tick label size
)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v1_e_max_diff.png", plot = pp_check_plots_m5v1_e_max_diff, width = 3, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v1_Filter_max_diff.png", plot = pp_check_plots_m5v1_Filter_max_diff, width = 3, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v1_e_prop_diff.png", plot = pp_check_plots_m5v1_e_prop_diff, width = 3, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v1_R_c_diff.png", plot = pp_check_plots_m5v1_R_c_diff, width = 3, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v1_G_c_diff.png", plot = pp_check_plots_m5v1_G_c_diff, width = 3, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v1_B_c_diff.png", plot = pp_check_plots_m5v1_B_c_diff, width = 3, height = 3, units = "in", dpi = 300)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
Warning message:
"The following arguments were unrecognized and ignored: nreps"
Coordinate system already present.
Adding new coordinate system, which
will replace the existing one.
Using 10 posterior draws for ppc type 'dens_overlay' by default.
Warning message:
"The following arguments were unrecognized and ignored: nreps"
Coordinate system already present.
Adding new coordinate system, which
will replace the existing one.
Using 10 posterior draws for ppc type 'dens_overlay' by default.
Warning message:
"The following arguments were unrecognized and ignored: nreps"
Coordinate system already present.
Adding new coordinate system, which
will replace the existing one.
Using 10 posterior draws for ppc type 'dens_overlay' by default.
Warning message:
"The following arguments were unrecognized and ignored: nreps"
Coordinate system already present.
Adding new coordinate system, which
will replace the existing one.
Using 10 posterior draws for ppc type 'dens_overlay' by default.
Warning message:
"The following arguments were unrecognized and ignored: nreps"
Coordinate system already present.
Adding new coordinate system, which
will replace the existing one.
Using 10 posterior draws for ppc type 'dens_overlay' by default.
Warning message:
"The following arguments were unrecognized and ignored: nreps"
Coordinate system already present.
Adding new coordinate system, which
will replace the existing one.
Show code cell source
# Convert images to base64
pp_check_plots_m5v1_e_max_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v1_e_max_diff.png")
pp_check_plots_m5v1_Filter_max_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v1_Filter_max_diff.png")
pp_check_plots_m5v1_e_prop_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v1_e_prop_diff.png")
pp_check_plots_m5v1_R_c_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v1_R_c_diff.png")
pp_check_plots_m5v1_G_c_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v1_G_c_diff.png")
pp_check_plots_m5v1_B_c_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v1_B_c_diff.png")
# Create the HTML (horizontal display)
html_pp_check_plots_m5v1 <- paste0("
<style>
.grid-container {
display: grid;
grid-template-columns: repeat(3, 1fr); /* 3 columns per row */
gap: 15px;
padding: 10px;
justify-items: center;
}
.grid-container img {
max-width: 300px;
width: 100%;
height: auto;
border: 1px solid #ccc;
}
</style>
<div class='grid-container'>
<img src='", pp_check_plots_m5v1_e_max_diff, "' alt='e_max Plot (Model 1)'>
<img src='", pp_check_plots_m5v1_Filter_max_diff, "' alt='Filter_max Plot (Model 1)'>
<img src='", pp_check_plots_m5v1_e_prop_diff, "' alt='e_prop Plot (Model 1)'>
<img src='", pp_check_plots_m5v1_R_c_diff, "' alt='R Reflectance Plot (Model 1)'>
<img src='", pp_check_plots_m5v1_G_c_diff, "' alt='G Reflectance Plot (Model 1)'>
<img src='", pp_check_plots_m5v1_B_c_diff, "' alt='B Reflectance Plot (Model 1)'>
</div>
")
IRdisplay::display_html(html_pp_check_plots_m5v1)
Show code cell content
# Generate posterior predictive check plots for all models
pp_check_plots_m5v2_e_max_diff <- generate_pp_check(
model = m5v2_e_max_diff,
nreps = 100,
axis_title_y = TRUE,
y_label = "Density",
plot_title = "e_max Background Matching Model",
axis_title_size = 10, # custom title size
axis_text_size = 8 # custom tick label size
)
pp_check_plots_m5v2_Filter_max_diff <- generate_pp_check(
model = m5v2_Filter_max_diff,
nreps = 100,
axis_title_y = TRUE,
y_label = "Density",
plot_title = "Filter_max Background Matching Model",
axis_title_size = 10, # custom title size
axis_text_size = 8 # custom tick label size
)
# Generate posterior predictive check plots for all models
pp_check_plots_m5v2_e_prop_diff <- generate_pp_check(
model = m5v2_e_prop_diff,
nreps = 100,
axis_title_y = TRUE,
y_label = "Density",
plot_title = "e_prop Background Matching Model",
axis_title_size = 10, # custom title size
axis_text_size = 8 # custom tick label size
)
pp_check_plots_m5v2_R_c_diff <- generate_pp_check(
model = m5v2_R_c_diff,
nreps = 100,
axis_title_y = TRUE,
y_label = "Density",
plot_title = "R Reflectance Background Matching Model",
axis_title_size = 10, # custom title size
axis_text_size = 8 # custom tick label size
)
# Generate posterior predictive check plots for all models
pp_check_plots_m5v2_G_c_diff <- generate_pp_check(
model = m5v2_G_c_diff,
nreps = 100,
axis_title_y = TRUE,
y_label = "Density",
plot_title = "G Reflectance Background Matching Model",
axis_title_size = 10, # custom title size
axis_text_size = 8 # custom tick label size
)
pp_check_plots_m5v2_B_c_diff <- generate_pp_check(
model = m5v2_B_c_diff,
nreps = 100,
axis_title_y = TRUE,
y_label = "Density",
plot_title = "B Reflectance Background Matching Model",
axis_title_size = 10, # custom title size
axis_text_size = 8 # custom tick label size
)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v2_e_max_diff.png", plot = pp_check_plots_m5v2_e_max_diff, width = 3, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v2_Filter_max_diff.png", plot = pp_check_plots_m5v2_Filter_max_diff, width = 3, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v2_e_prop_diff.png", plot = pp_check_plots_m5v2_e_prop_diff, width = 3, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v2_R_c_diff.png", plot = pp_check_plots_m5v2_R_c_diff, width = 3, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v2_G_c_diff.png", plot = pp_check_plots_m5v2_G_c_diff, width = 3, height = 3, units = "in", dpi = 300)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v2_B_c_diff.png", plot = pp_check_plots_m5v2_B_c_diff, width = 3, height = 3, units = "in", dpi = 300)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
Warning message:
"The following arguments were unrecognized and ignored: nreps"
Coordinate system already present.
Adding new coordinate system, which
will replace the existing one.
Using 10 posterior draws for ppc type 'dens_overlay' by default.
Warning message:
"The following arguments were unrecognized and ignored: nreps"
Coordinate system already present.
Adding new coordinate system, which
will replace the existing one.
Using 10 posterior draws for ppc type 'dens_overlay' by default.
Warning message:
"The following arguments were unrecognized and ignored: nreps"
Coordinate system already present.
Adding new coordinate system, which
will replace the existing one.
Using 10 posterior draws for ppc type 'dens_overlay' by default.
Warning message:
"The following arguments were unrecognized and ignored: nreps"
Coordinate system already present.
Adding new coordinate system, which
will replace the existing one.
Using 10 posterior draws for ppc type 'dens_overlay' by default.
Warning message:
"The following arguments were unrecognized and ignored: nreps"
Coordinate system already present.
Adding new coordinate system, which
will replace the existing one.
Using 10 posterior draws for ppc type 'dens_overlay' by default.
Warning message:
"The following arguments were unrecognized and ignored: nreps"
Coordinate system already present.
Adding new coordinate system, which
will replace the existing one.
Show code cell source
# Convert images to base64
pp_check_plots_m5v2_e_max_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v2_e_max_diff.png")
pp_check_plots_m5v2_Filter_max_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v2_Filter_max_diff.png")
pp_check_plots_m5v2_e_prop_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v2_e_prop_diff.png")
pp_check_plots_m5v2_R_c_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v2_R_c_diff.png")
pp_check_plots_m5v2_G_c_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v2_G_c_diff.png")
pp_check_plots_m5v2_B_c_diff <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/pp_check_plots_m5v2_B_c_diff.png")
# Create the HTML (horizontal display)
html_pp_check_plots_m5v2 <- paste0("
<style>
.grid-container {
display: grid;
grid-template-columns: repeat(3, 1fr); /* 3 columns per row */
gap: 15px;
padding: 10px;
justify-items: center;
}
.grid-container img {
max-width: 300px;
width: 100%;
height: auto;
border: 1px solid #ccc;
}
</style>
<div class='grid-container'>
<img src='", pp_check_plots_m5v2_e_max_diff, "' alt='e_max Plot (Model 2)'>
<img src='", pp_check_plots_m5v2_Filter_max_diff, "' alt='Filter_max Plot (Model 2)'>
<img src='", pp_check_plots_m5v2_e_prop_diff, "' alt='e_prop Plot (Model 2)'>
<img src='", pp_check_plots_m5v2_R_c_diff, "' alt='R Reflectance Plot (Model 2)'>
<img src='", pp_check_plots_m5v2_G_c_diff, "' alt='G Reflectance Plot (Model 2)'>
<img src='", pp_check_plots_m5v2_B_c_diff, "' alt='B Reflectance Plot (Model 2)'>
</div>
")
IRdisplay::display_html(html_pp_check_plots_m5v2)
Check convergence#
Check that all \(\hat{R}\) values are close to 1, indicating good convergence.
Show code cell content
# Create a helper function
extract_rhat <- function(model, model_name) {
rhat(model) %>%
as.data.frame() %>%
rownames_to_column(var = "Parameter") %>%
dplyr::filter(startsWith(Parameter, "b_")) %>% # <-- keep only b_ terms
dplyr::rename(Rhat = 2) %>%
dplyr::mutate(Model = model_name) %>%
dplyr::mutate(across(where(is.numeric), ~ signif(.x, digits = 3)))
}
# Extract for each model group
rhat1_m5v1 <- extract_rhat(m5v1_e_max_diff, "e_max_diff")
rhat2_m5v1 <- extract_rhat(m5v1_Filter_max_diff, "Filter_max_diff")
rhat3_m5v1 <- extract_rhat(m5v1_e_prop_diff, "e_prop_diff")
rhat4_m5v1 <- extract_rhat(m5v1_R_c_diff, "R_c_diff")
rhat5_m5v1 <- extract_rhat(m5v1_G_c_diff, "G_c_diff")
rhat6_m5v1 <- extract_rhat(m5v1_B_c_diff, "B_c_diff")
Show code cell content
# Save tables
write.table(rhat1_m5v1, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_rhat_m5v1_e_max_diff.csv", sep = ",", row.names = FALSE, col.names = TRUE)
write.table(rhat2_m5v1, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_rhat_m5v1_Filter_max_diff.csv", sep = ",", row.names = FALSE, col.names = TRUE)
write.table(rhat3_m5v1, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_rhat_m5v1_e_prop_diff.csv", sep = ",", row.names = FALSE, col.names = TRUE)
write.table(rhat4_m5v1, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_rhat_m5v1_R_c_diff.csv", sep = ",", row.names = FALSE, col.names = TRUE)
write.table(rhat5_m5v1, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_rhat_m5v1_G_c_diff.csv", sep = ",", row.names = FALSE, col.names = TRUE)
write.table(rhat6_m5v1, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_rhat_m5v1_B_c_diff.csv", sep = ",", row.names = FALSE, col.names = TRUE)
Show code cell source
# Convert each data frame to a plain HTML table string
table_1_rhat_m5v1 <- minimal_html_table(rhat1_m5v1, caption = "Rhat values - e_max Background Matching")
table_2_rhat_m5v1 <- minimal_html_table(rhat2_m5v1, caption = "Rhat values - Filter_max Background Matching")
table_3_rhat_m5v1 <- minimal_html_table(rhat3_m5v1, caption = "Rhat values - e_prop Background Matching")
table_4_rhat_m5v1 <- minimal_html_table(rhat4_m5v1, caption = "Rhat values - R Background Matching")
table_5_rhat_m5v1 <- minimal_html_table(rhat5_m5v1, caption = "Rhat values - G Background Matching")
table_6_rhat_m5v1 <- minimal_html_table(rhat6_m5v1, caption = "Rhat values - B Background Matching")
my_tabs_rhat_m5v1 <- '
<style>
/* Basic container styling */
.tabs-container {
width: 100%;
margin: 1em 0;
}
/* Hide the radio inputs (we only show their labels as tabs) */
.tabs-container input[type="radio"] {
display: none;
}
/* The “tab-label” styling: looks like a tab */
.tab-label {
display: inline-block;
padding: 10px;
margin-right: 2px;
background: #eee;
border: 1px solid #ccc;
cursor: pointer;
border-bottom: none;
}
/* The active tab label */
.tab-label-active {
background: #fff;
}
/* The panel that holds table content */
.tab-content {
border: 1px solid #ccc;
padding: 10px;
display: none;
}
/* For each radio input, show its corresponding content when checked */
#tab1_rhat_m5v1:checked ~ #content1_rhat_m5v1,
#tab2_rhat_m5v1:checked ~ #content2_rhat_m5v1,
#tab3_rhat_m5v1:checked ~ #content3_rhat_m5v1,
#tab4_rhat_m5v1:checked ~ #content4_rhat_m5v1,
#tab5_rhat_m5v1:checked ~ #content5_rhat_m5v1,
#tab6_rhat_m5v1:checked ~ #content6_rhat_m5v1 {
display: block;
}
/* Also style the label of the checked radio as “active” using the :checked + label technique */
#tab1_rhat_m5v1:checked + label[for="tab1_rhat_m5v1"],
#tab2_rhat_m5v1:checked + label[for="tab2_rhat_m5v1"],
#tab3_rhat_m5v1:checked + label[for="tab3_rhat_m5v1"],
#tab4_rhat_m5v1:checked + label[for="tab4_rhat_m5v1"],
#tab5_rhat_m5v1:checked + label[for="tab5_rhat_m5v1"],
#tab6_rhat_m5v1:checked + label[for="tab6_rhat_m5v1"] {
background: #fff;
border-bottom: none;
}
</style>
<div class="tabs-container">
<!-- 1) Tab radio + label -->
<input type="radio" name="tabs_rhat_m5v1" id="tab1_rhat_m5v1" checked>
<label class="tab-label" for="tab1_rhat_m5v1">Table 1</label>
<!-- 2) Tab radio + label -->
<input type="radio" name="tabs_rhat_m5v1" id="tab2_rhat_m5v1">
<label class="tab-label" for="tab2_rhat_m5v1">Table 2</label>
<!-- 3) Tab radio + label -->
<input type="radio" name="tabs_rhat_m5v1" id="tab3_rhat_m5v1">
<label class="tab-label" for="tab3_rhat_m5v1">Table 3</label>
<!-- 4) Tab radio + label -->
<input type="radio" name="tabs_rhat_m5v1" id="tab4_rhat_m5v1">
<label class="tab-label" for="tab4_rhat_m5v1">Table 4</label>
<!-- 5) Tab radio + label -->
<input type="radio" name="tabs_rhat_m5v1" id="tab5_rhat_m5v1">
<label class="tab-label" for="tab5_rhat_m5v1">Table 5</label>
<!-- 6) Tab radio + label -->
<input type="radio" name="tabs_rhat_m5v1" id="tab6_rhat_m5v1">
<label class="tab-label" for="tab6_rhat_m5v1">Table 6</label>
<!-- Content for each tab -->
<div class="tab-content" id="content1_rhat_m5v1">REPLACE_WITH_table_1_m5v1</div>
<div class="tab-content" id="content2_rhat_m5v1">REPLACE_WITH_table_2_m5v1</div>
<div class="tab-content" id="content3_rhat_m5v1">REPLACE_WITH_table_3_m5v1</div>
<div class="tab-content" id="content4_rhat_m5v1">REPLACE_WITH_table_4_m5v1</div>
<div class="tab-content" id="content5_rhat_m5v1">REPLACE_WITH_table_5_m5v1</div>
<div class="tab-content" id="content6_rhat_m5v1">REPLACE_WITH_table_6_m5v1</div>
</div>
'
# Now do the replacements for each table
my_tabs_rhat_m5v1 <- gsub("REPLACE_WITH_table_1_m5v1", table_1_rhat_m5v1, my_tabs_rhat_m5v1)
my_tabs_rhat_m5v1 <- gsub("REPLACE_WITH_table_2_m5v1", table_2_rhat_m5v1, my_tabs_rhat_m5v1)
my_tabs_rhat_m5v1 <- gsub("REPLACE_WITH_table_3_m5v1", table_3_rhat_m5v1, my_tabs_rhat_m5v1)
my_tabs_rhat_m5v1 <- gsub("REPLACE_WITH_table_4_m5v1", table_4_rhat_m5v1, my_tabs_rhat_m5v1)
my_tabs_rhat_m5v1 <- gsub("REPLACE_WITH_table_5_m5v1", table_5_rhat_m5v1, my_tabs_rhat_m5v1)
my_tabs_rhat_m5v1 <- gsub("REPLACE_WITH_table_6_m5v1", table_6_rhat_m5v1, my_tabs_rhat_m5v1)
IRdisplay::display_html(my_tabs_rhat_m5v1)
| Parameter | Rhat | Model |
|---|---|---|
| b_Intercept | 1 | e_max_diff |
| b_SexF | 1 | e_max_diff |
| b_MicrohabitatHydroid_Bryozoa | 1 | e_max_diff |
| b_MicrohabitatRed_Algae | 1 | e_max_diff |
| b_Viewpoint10 | 1 | e_max_diff |
| b_Viewpoint20 | 1 | e_max_diff |
| Parameter | Rhat | Model |
|---|---|---|
| b_Intercept | 1 | Filter_max_diff |
| b_SexF | 1 | Filter_max_diff |
| b_MicrohabitatHydroid_Bryozoa | 1 | Filter_max_diff |
| b_MicrohabitatRed_Algae | 1 | Filter_max_diff |
| b_Viewpoint10 | 1 | Filter_max_diff |
| b_Viewpoint20 | 1 | Filter_max_diff |
| Parameter | Rhat | Model |
|---|---|---|
| b_Intercept | 1 | e_prop_diff |
| b_SexF | 1 | e_prop_diff |
| b_MicrohabitatHydroid_Bryozoa | 1 | e_prop_diff |
| b_MicrohabitatRed_Algae | 1 | e_prop_diff |
| b_Viewpoint10 | 1 | e_prop_diff |
| b_Viewpoint20 | 1 | e_prop_diff |
| Parameter | Rhat | Model |
|---|---|---|
| b_Intercept | 1.01 | R_c_diff |
| b_SexF | 1 | R_c_diff |
| b_MicrohabitatHydroid_Bryozoa | 1 | R_c_diff |
| b_MicrohabitatRed_Algae | 1 | R_c_diff |
| b_Viewpoint10 | 1 | R_c_diff |
| b_Viewpoint20 | 1 | R_c_diff |
| Parameter | Rhat | Model |
|---|---|---|
| b_Intercept | 1 | G_c_diff |
| b_SexF | 1 | G_c_diff |
| b_MicrohabitatHydroid_Bryozoa | 1 | G_c_diff |
| b_MicrohabitatRed_Algae | 1 | G_c_diff |
| b_Viewpoint10 | 1 | G_c_diff |
| b_Viewpoint20 | 1 | G_c_diff |
| Parameter | Rhat | Model |
|---|---|---|
| b_Intercept | 1 | B_c_diff |
| b_SexF | 1 | B_c_diff |
| b_MicrohabitatHydroid_Bryozoa | 1 | B_c_diff |
| b_MicrohabitatRed_Algae | 1 | B_c_diff |
| b_Viewpoint10 | 1 | B_c_diff |
| b_Viewpoint20 | 1 | B_c_diff |
Show code cell content
# Create a helper function
extract_rhat <- function(model, model_name) {
rhat(model) %>%
as.data.frame() %>%
rownames_to_column(var = "Parameter") %>%
dplyr::filter(startsWith(Parameter, "b_")) %>% # <-- keep only b_ terms
dplyr::rename(Rhat = 2) %>%
dplyr::mutate(Model = model_name) %>%
dplyr::mutate(across(where(is.numeric), ~ signif(.x, digits = 3)))
}
# Extract for each model group
rhat1_m5v2 <- extract_rhat(m5v2_e_max_diff, "e_max_diff")
rhat2_m5v2 <- extract_rhat(m5v2_Filter_max_diff, "Filter_max_diff")
rhat3_m5v2 <- extract_rhat(m5v2_e_prop_diff, "e_prop_diff")
rhat4_m5v2 <- extract_rhat(m5v2_R_c_diff, "R_c_diff")
rhat5_m5v2 <- extract_rhat(m5v2_G_c_diff, "G_c_diff")
rhat6_m5v2 <- extract_rhat(m5v2_B_c_diff, "B_c_diff")
Show code cell content
# Save tables
write.table(rhat1_m5v2, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_rhat_m5v2_e_max_diff.csv", sep = ",", row.names = FALSE, col.names = TRUE)
write.table(rhat2_m5v2, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_rhat_m5v2_Filter_max_diff.csv", sep = ",", row.names = FALSE, col.names = TRUE)
write.table(rhat3_m5v2, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_rhat_m5v2_e_prop_diff.csv", sep = ",", row.names = FALSE, col.names = TRUE)
write.table(rhat4_m5v2, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_rhat_m5v2_R_c_diff.csv", sep = ",", row.names = FALSE, col.names = TRUE)
write.table(rhat5_m5v2, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_rhat_m5v2_G_c_diff.csv", sep = ",", row.names = FALSE, col.names = TRUE)
write.table(rhat6_m5v2, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_rhat_m5v2_B_c_diff.csv", sep = ",", row.names = FALSE, col.names = TRUE)
Show code cell source
# Convert each data frame to a plain HTML table string
table_1_rhat_m5v2 <- minimal_html_table(rhat1_m5v2, caption = "Rhat values - e_max Background Matching")
table_2_rhat_m5v2 <- minimal_html_table(rhat2_m5v2, caption = "Rhat values - Filter_max Background Matching")
table_3_rhat_m5v2 <- minimal_html_table(rhat3_m5v2, caption = "Rhat values - e_prop Background Matching")
table_4_rhat_m5v2 <- minimal_html_table(rhat4_m5v2, caption = "Rhat values - R Background Matching")
table_5_rhat_m5v2 <- minimal_html_table(rhat5_m5v2, caption = "Rhat values - G Background Matching")
table_6_rhat_m5v2 <- minimal_html_table(rhat6_m5v2, caption = "Rhat values - B Background Matching")
my_tabs_rhat_m5v2 <- '
<style>
/* Basic container styling */
.tabs-container {
width: 100%;
margin: 1em 0;
}
/* Hide the radio inputs (we only show their labels as tabs) */
.tabs-container input[type="radio"] {
display: none;
}
/* The “tab-label” styling: looks like a tab */
.tab-label {
display: inline-block;
padding: 10px;
margin-right: 2px;
background: #eee;
border: 1px solid #ccc;
cursor: pointer;
border-bottom: none;
}
/* The active tab label */
.tab-label-active {
background: #fff;
}
/* The panel that holds table content */
.tab-content {
border: 1px solid #ccc;
padding: 10px;
display: none;
}
/* For each radio input, show its corresponding content when checked */
#tab1_rhat_m5v2:checked ~ #content1_rhat_m5v2,
#tab2_rhat_m5v2:checked ~ #content2_rhat_m5v2,
#tab3_rhat_m5v2:checked ~ #content3_rhat_m5v2,
#tab4_rhat_m5v2:checked ~ #content4_rhat_m5v2,
#tab5_rhat_m5v2:checked ~ #content5_rhat_m5v2,
#tab6_rhat_m5v2:checked ~ #content6_rhat_m5v2 {
display: block;
}
/* Also style the label of the checked radio as “active” using the :checked + label technique */
#tab1_rhat_m5v2:checked + label[for="tab1_rhat_m5v2"],
#tab2_rhat_m5v2:checked + label[for="tab2_rhat_m5v2"],
#tab3_rhat_m5v2:checked + label[for="tab3_rhat_m5v2"],
#tab4_rhat_m5v2:checked + label[for="tab4_rhat_m5v2"],
#tab5_rhat_m5v2:checked + label[for="tab5_rhat_m5v2"],
#tab6_rhat_m5v2:checked + label[for="tab6_rhat_m5v2"] {
background: #fff;
border-bottom: none;
}
</style>
<div class="tabs-container">
<!-- 1) Tab radio + label -->
<input type="radio" name="tabs_rhat_m5v2" id="tab1_rhat_m5v2" checked>
<label class="tab-label" for="tab1_rhat_m5v2">Table 1</label>
<!-- 2) Tab radio + label -->
<input type="radio" name="tabs_rhat_m5v2" id="tab2_rhat_m5v2">
<label class="tab-label" for="tab2_rhat_m5v2">Table 2</label>
<!-- 3) Tab radio + label -->
<input type="radio" name="tabs_rhat_m5v2" id="tab3_rhat_m5v2">
<label class="tab-label" for="tab3_rhat_m5v2">Table 3</label>
<!-- 4) Tab radio + label -->
<input type="radio" name="tabs_rhat_m5v2" id="tab4_rhat_m5v2">
<label class="tab-label" for="tab4_rhat_m5v2">Table 4</label>
<!-- 5) Tab radio + label -->
<input type="radio" name="tabs_rhat_m5v2" id="tab5_rhat_m5v2">
<label class="tab-label" for="tab5_rhat_m5v2">Table 5</label>
<!-- 6) Tab radio + label -->
<input type="radio" name="tabs_rhat_m5v2" id="tab6_rhat_m5v2">
<label class="tab-label" for="tab6_rhat_m5v2">Table 6</label>
<!-- Content for each tab -->
<div class="tab-content" id="content1_rhat_m5v2">REPLACE_WITH_table_1_m5v2</div>
<div class="tab-content" id="content2_rhat_m5v2">REPLACE_WITH_table_2_m5v2</div>
<div class="tab-content" id="content3_rhat_m5v2">REPLACE_WITH_table_3_m5v2</div>
<div class="tab-content" id="content4_rhat_m5v2">REPLACE_WITH_table_4_m5v2</div>
<div class="tab-content" id="content5_rhat_m5v2">REPLACE_WITH_table_5_m5v2</div>
<div class="tab-content" id="content6_rhat_m5v2">REPLACE_WITH_table_6_m5v2</div>
</div>
'
# Now do the replacements for each table
my_tabs_rhat_m5v2 <- gsub("REPLACE_WITH_table_1_m5v2", table_1_rhat_m5v2, my_tabs_rhat_m5v2)
my_tabs_rhat_m5v2 <- gsub("REPLACE_WITH_table_2_m5v2", table_2_rhat_m5v2, my_tabs_rhat_m5v2)
my_tabs_rhat_m5v2 <- gsub("REPLACE_WITH_table_3_m5v2", table_3_rhat_m5v2, my_tabs_rhat_m5v2)
my_tabs_rhat_m5v2 <- gsub("REPLACE_WITH_table_4_m5v2", table_4_rhat_m5v2, my_tabs_rhat_m5v2)
my_tabs_rhat_m5v2 <- gsub("REPLACE_WITH_table_5_m5v2", table_5_rhat_m5v2, my_tabs_rhat_m5v2)
my_tabs_rhat_m5v2 <- gsub("REPLACE_WITH_table_6_m5v2", table_6_rhat_m5v2, my_tabs_rhat_m5v2)
IRdisplay::display_html(my_tabs_rhat_m5v2)
| Parameter | Rhat | Model |
|---|---|---|
| b_Intercept | 1 | e_max_diff |
| b_SexF | 1 | e_max_diff |
| b_Microhabitat_AssociationPresent | 1 | e_max_diff |
| b_Viewpoint10 | 1 | e_max_diff |
| b_Viewpoint20 | 1 | e_max_diff |
| Parameter | Rhat | Model |
|---|---|---|
| b_Intercept | 1 | Filter_max_diff |
| b_SexF | 1 | Filter_max_diff |
| b_Microhabitat_AssociationPresent | 1 | Filter_max_diff |
| b_Viewpoint10 | 1 | Filter_max_diff |
| b_Viewpoint20 | 1 | Filter_max_diff |
| Parameter | Rhat | Model |
|---|---|---|
| b_Intercept | 1 | e_prop_diff |
| b_SexF | 1 | e_prop_diff |
| b_Microhabitat_AssociationPresent | 1 | e_prop_diff |
| b_Viewpoint10 | 1 | e_prop_diff |
| b_Viewpoint20 | 1 | e_prop_diff |
| Parameter | Rhat | Model |
|---|---|---|
| b_Intercept | 1 | R_c_diff |
| b_SexF | 1 | R_c_diff |
| b_Microhabitat_AssociationPresent | 1 | R_c_diff |
| b_Viewpoint10 | 1 | R_c_diff |
| b_Viewpoint20 | 1 | R_c_diff |
| Parameter | Rhat | Model |
|---|---|---|
| b_Intercept | 1.01 | G_c_diff |
| b_SexF | 1 | G_c_diff |
| b_Microhabitat_AssociationPresent | 1 | G_c_diff |
| b_Viewpoint10 | 1 | G_c_diff |
| b_Viewpoint20 | 1 | G_c_diff |
| Parameter | Rhat | Model |
|---|---|---|
| b_Intercept | 1 | B_c_diff |
| b_SexF | 1 | B_c_diff |
| b_Microhabitat_AssociationPresent | 1 | B_c_diff |
| b_Viewpoint10 | 1 | B_c_diff |
| b_Viewpoint20 | 1 | B_c_diff |
Check uncertainty#
Extract parameter estimates and their confidence intervals to assess the significance of the predictors on color pattern metrics. We check 85% and 95% confidence intervals. Summaries are displayed in tables for all models.
Show code cell content
# Extract summaries for each variable
extract_summary <- function(model, prob_85, prob_95) {
summary_85 <- summary(model, prob = prob_85)
summary_95 <- summary(model, prob = prob_95)
as.data.frame(summary_85$fixed) %>%
dplyr::select("Estimate", "Est.Error", "l-85% CI", "u-85% CI") %>%
mutate(
"l-95% CI" = summary_95$fixed$`l-95% CI`,
"u-95% CI" = summary_95$fixed$`u-95% CI`
) %>%
mutate(across(where(is.numeric), ~ signif(.x, digits = 3))) %>%
rownames_to_column(var = "Parameter") # Add rownames as Parameter column
}
uncertainty1_m5v1 <- extract_summary(m5v1_e_max_diff, prob_85 = 0.85, prob_95 = 0.95)
uncertainty2_m5v1 <- extract_summary(m5v1_Filter_max_diff, prob_85 = 0.85, prob_95 = 0.95)
uncertainty3_m5v1 <- extract_summary(m5v1_e_prop_diff, prob_85 = 0.85, prob_95 = 0.95)
uncertainty4_m5v1 <- extract_summary(m5v1_R_c_diff, prob_85 = 0.85, prob_95 = 0.95)
uncertainty5_m5v1 <- extract_summary(m5v1_G_c_diff, prob_85 = 0.85, prob_95 = 0.95)
uncertainty6_m5v1 <- extract_summary(m5v1_B_c_diff, prob_85 = 0.85, prob_95 = 0.95)
Show code cell content
# Save tables
write.table(uncertainty1_m5v1, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_uncertainty_m5v1_e_max_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(uncertainty2_m5v1, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_uncertainty_m5v1_Filter_max_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(uncertainty3_m5v1, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_uncertainty_m5v1_e_prop_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(uncertainty4_m5v1, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_uncertainty_m5v1_R_c_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(uncertainty5_m5v1, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_uncertainty_m5v1_G_c_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(uncertainty6_m5v1, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_uncertainty_m5v1_B_c_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
Show code cell source
# Convert each data frame to a plain HTML table string
table_1_uncertainty_m5v1 <- minimal_html_table(uncertainty1_m5v1, caption = "Uncertainty values - e_max Background Matching")
table_2_uncertainty_m5v1 <- minimal_html_table(uncertainty2_m5v1, caption = "Uncertainty values - Filter_max Background Matching")
table_3_uncertainty_m5v1 <- minimal_html_table(uncertainty3_m5v1, caption = "Uncertainty values - e_prop Background Matching")
table_4_uncertainty_m5v1 <- minimal_html_table(uncertainty4_m5v1, caption = "Uncertainty values - R Background Matching")
table_5_uncertainty_m5v1 <- minimal_html_table(uncertainty5_m5v1, caption = "Uncertainty values - G Background Matching")
table_6_uncertainty_m5v1 <- minimal_html_table(uncertainty6_m5v1, caption = "Uncertainty values - B Background Matching")
my_tabs_uncertainty_m5v1 <- '
<style>
/* Basic container styling */
.tabs-container {
width: 100%;
margin: 1em 0;
}
/* Hide the radio inputs (we only show their labels as tabs) */
.tabs-container input[type="radio"] {
display: none;
}
/* The “tab-label” styling: looks like a tab */
.tab-label {
display: inline-block;
padding: 10px;
margin-right: 2px;
background: #eee;
border: 1px solid #ccc;
cursor: pointer;
border-bottom: none;
}
/* The active tab label */
.tab-label-active {
background: #fff;
}
/* The panel that holds table content */
.tab-content {
border: 1px solid #ccc;
padding: 10px;
display: none;
}
/* For each radio input, show its corresponding content when checked */
#tab1_uncertainty_m5v1:checked ~ #content1_uncertainty_m5v1,
#tab2_uncertainty_m5v1:checked ~ #content2_uncertainty_m5v1,
#tab3_uncertainty_m5v1:checked ~ #content3_uncertainty_m5v1,
#tab4_uncertainty_m5v1:checked ~ #content4_uncertainty_m5v1,
#tab5_uncertainty_m5v1:checked ~ #content5_uncertainty_m5v1,
#tab6_uncertainty_m5v1:checked ~ #content6_uncertainty_m5v1 {
display: block;
}
/* Also style the label of the checked radio as “active” using the :checked + label technique */
#tab1_uncertainty_m5v1:checked + label[for="tab1_uncertainty_m5v1"],
#tab2_uncertainty_m5v1:checked + label[for="tab2_uncertainty_m5v1"],
#tab3_uncertainty_m5v1:checked + label[for="tab3_uncertainty_m5v1"],
#tab4_uncertainty_m5v1:checked + label[for="tab4_uncertainty_m5v1"],
#tab5_uncertainty_m5v1:checked + label[for="tab5_uncertainty_m5v1"],
#tab6_uncertainty_m5v1:checked + label[for="tab6_uncertainty_m5v1"] {
background: #fff;
border-bottom: none;
}
</style>
<div class="tabs-container">
<!-- 1) Tab radio + label -->
<input type="radio" name="tabs_uncertainty_m5v1" id="tab1_uncertainty_m5v1" checked>
<label class="tab-label" for="tab1_uncertainty_m5v1">Table 1</label>
<!-- 2) Tab radio + label -->
<input type="radio" name="tabs_uncertainty_m5v1" id="tab2_uncertainty_m5v1">
<label class="tab-label" for="tab2_uncertainty_m5v1">Table 2</label>
<!-- 3) Tab radio + label -->
<input type="radio" name="tabs_uncertainty_m5v1" id="tab3_uncertainty_m5v1">
<label class="tab-label" for="tab3_uncertainty_m5v1">Table 3</label>
<!-- 4) Tab radio + label -->
<input type="radio" name="tabs_uncertainty_m5v1" id="tab4_uncertainty_m5v1">
<label class="tab-label" for="tab4_uncertainty_m5v1">Table 4</label>
<!-- 5) Tab radio + label -->
<input type="radio" name="tabs_uncertainty_m5v1" id="tab5_uncertainty_m5v1">
<label class="tab-label" for="tab5_uncertainty_m5v1">Table 5</label>
<!-- 6) Tab radio + label -->
<input type="radio" name="tabs_uncertainty_m5v1" id="tab6_uncertainty_m5v1">
<label class="tab-label" for="tab6_uncertainty_m5v1">Table 6</label>
<!-- Content for each tab -->
<div class="tab-content" id="content1_uncertainty_m5v1">REPLACE_WITH_table_1_m5v1</div>
<div class="tab-content" id="content2_uncertainty_m5v1">REPLACE_WITH_table_2_m5v1</div>
<div class="tab-content" id="content3_uncertainty_m5v1">REPLACE_WITH_table_3_m5v1</div>
<div class="tab-content" id="content4_uncertainty_m5v1">REPLACE_WITH_table_4_m5v1</div>
<div class="tab-content" id="content5_uncertainty_m5v1">REPLACE_WITH_table_5_m5v1</div>
<div class="tab-content" id="content6_uncertainty_m5v1">REPLACE_WITH_table_6_m5v1</div>
</div>
'
# Now do the replacements for each table
my_tabs_uncertainty_m5v1 <- gsub("REPLACE_WITH_table_1_m5v1", table_1_uncertainty_m5v1, my_tabs_uncertainty_m5v1)
my_tabs_uncertainty_m5v1 <- gsub("REPLACE_WITH_table_2_m5v1", table_2_uncertainty_m5v1, my_tabs_uncertainty_m5v1)
my_tabs_uncertainty_m5v1 <- gsub("REPLACE_WITH_table_3_m5v1", table_3_uncertainty_m5v1, my_tabs_uncertainty_m5v1)
my_tabs_uncertainty_m5v1 <- gsub("REPLACE_WITH_table_4_m5v1", table_4_uncertainty_m5v1, my_tabs_uncertainty_m5v1)
my_tabs_uncertainty_m5v1 <- gsub("REPLACE_WITH_table_5_m5v1", table_5_uncertainty_m5v1, my_tabs_uncertainty_m5v1)
my_tabs_uncertainty_m5v1 <- gsub("REPLACE_WITH_table_6_m5v1", table_6_uncertainty_m5v1, my_tabs_uncertainty_m5v1)
IRdisplay::display_html(my_tabs_uncertainty_m5v1)
| Parameter | Estimate | Est.Error | l-85% CI | u-85% CI | l-95% CI | u-95% CI |
|---|---|---|---|---|---|---|
| Intercept | -0.776 | 0.0334 | -0.824 | -0.729 | -0.841 | -0.71 |
| SexF | 0.105 | 0.0117 | 0.0884 | 0.122 | 0.0824 | 0.128 |
| MicrohabitatHydroid_Bryozoa | 0.0534 | 0.0134 | 0.0343 | 0.0726 | 0.0272 | 0.0795 |
| MicrohabitatRed_Algae | -0.198 | 0.0187 | -0.225 | -0.171 | -0.235 | -0.162 |
| Viewpoint10 | -0.0424 | 0.0136 | -0.062 | -0.0228 | -0.0693 | -0.0159 |
| Viewpoint20 | -0.0668 | 0.0138 | -0.0866 | -0.0469 | -0.0938 | -0.0399 |
| Parameter | Estimate | Est.Error | l-85% CI | u-85% CI | l-95% CI | u-95% CI |
|---|---|---|---|---|---|---|
| Intercept | 2.64 | 0.0198 | 2.61 | 2.66 | 2.6 | 2.67 |
| SexF | 0.0395 | 0.00766 | 0.0284 | 0.0506 | 0.0245 | 0.0548 |
| MicrohabitatHydroid_Bryozoa | -0.0529 | 0.00863 | -0.0653 | -0.0406 | -0.0699 | -0.0361 |
| MicrohabitatRed_Algae | 0.0578 | 0.0123 | 0.04 | 0.0755 | 0.0335 | 0.0825 |
| Viewpoint10 | -1.03 | 0.00907 | -1.04 | -1.02 | -1.05 | -1.01 |
| Viewpoint20 | -1.32 | 0.00907 | -1.33 | -1.31 | -1.34 | -1.3 |
| Parameter | Estimate | Est.Error | l-85% CI | u-85% CI | l-95% CI | u-95% CI |
|---|---|---|---|---|---|---|
| Intercept | -4.07 | 0.0201 | -4.1 | -4.05 | -4.11 | -4.03 |
| SexF | -0.246 | 0.00885 | -0.259 | -0.234 | -0.264 | -0.229 |
| MicrohabitatHydroid_Bryozoa | -0.0338 | 0.01 | -0.0482 | -0.0194 | -0.0534 | -0.0141 |
| MicrohabitatRed_Algae | 0.262 | 0.0142 | 0.241 | 0.282 | 0.234 | 0.289 |
| Viewpoint10 | 0.0817 | 0.0106 | 0.0662 | 0.0969 | 0.0607 | 0.103 |
| Viewpoint20 | 0.227 | 0.0106 | 0.212 | 0.242 | 0.206 | 0.248 |
| Parameter | Estimate | Est.Error | l-85% CI | u-85% CI | l-95% CI | u-95% CI |
|---|---|---|---|---|---|---|
| Intercept | 1.58 | 0.0411 | 1.52 | 1.64 | 1.5 | 1.67 |
| SexF | 0.175 | 0.0107 | 0.16 | 0.19 | 0.154 | 0.196 |
| MicrohabitatHydroid_Bryozoa | 0.0751 | 0.0121 | 0.0578 | 0.0925 | 0.0517 | 0.0991 |
| MicrohabitatRed_Algae | -0.0973 | 0.0172 | -0.122 | -0.0723 | -0.131 | -0.0633 |
| Viewpoint10 | 0.00615 | 0.0126 | -0.012 | 0.0244 | -0.0187 | 0.0307 |
| Viewpoint20 | 0.00782 | 0.0126 | -0.0106 | 0.0259 | -0.0171 | 0.0321 |
| Parameter | Estimate | Est.Error | l-85% CI | u-85% CI | l-95% CI | u-95% CI |
|---|---|---|---|---|---|---|
| Intercept | 1.53 | 0.0361 | 1.48 | 1.58 | 1.46 | 1.6 |
| SexF | 0.149 | 0.0114 | 0.132 | 0.165 | 0.126 | 0.17 |
| MicrohabitatHydroid_Bryozoa | 0.113 | 0.013 | 0.0942 | 0.131 | 0.0872 | 0.139 |
| MicrohabitatRed_Algae | -0.159 | 0.0184 | -0.186 | -0.133 | -0.195 | -0.123 |
| Viewpoint10 | 0.006 | 0.0134 | -0.0133 | 0.0253 | -0.0203 | 0.0325 |
| Viewpoint20 | 0.00661 | 0.0134 | -0.0127 | 0.0259 | -0.0197 | 0.0329 |
| Parameter | Estimate | Est.Error | l-85% CI | u-85% CI | l-95% CI | u-95% CI |
|---|---|---|---|---|---|---|
| Intercept | 1.22 | 0.0382 | 1.17 | 1.28 | 1.15 | 1.3 |
| SexF | 0.243 | 0.0129 | 0.224 | 0.261 | 0.218 | 0.268 |
| MicrohabitatHydroid_Bryozoa | 0.109 | 0.0147 | 0.0881 | 0.13 | 0.0806 | 0.138 |
| MicrohabitatRed_Algae | -0.34 | 0.0211 | -0.37 | -0.309 | -0.381 | -0.298 |
| Viewpoint10 | 0.00624 | 0.0154 | -0.0159 | 0.0284 | -0.0239 | 0.0363 |
| Viewpoint20 | 0.00677 | 0.0154 | -0.0154 | 0.029 | -0.0236 | 0.0369 |
Show code cell content
# Extract summaries for each variable
extract_summary <- function(model, prob_85, prob_95) {
summary_85 <- summary(model, prob = prob_85)
summary_95 <- summary(model, prob = prob_95)
as.data.frame(summary_85$fixed) %>%
dplyr::select("Estimate", "Est.Error", "l-85% CI", "u-85% CI") %>%
mutate(
"l-95% CI" = summary_95$fixed$`l-95% CI`,
"u-95% CI" = summary_95$fixed$`u-95% CI`
) %>%
mutate(across(where(is.numeric), ~ signif(.x, digits = 3))) %>%
rownames_to_column(var = "Parameter") # Add rownames as Parameter column
}
uncertainty1_m5v2 <- extract_summary(m5v2_e_max_diff, prob_85 = 0.85, prob_95 = 0.95)
uncertainty2_m5v2 <- extract_summary(m5v2_Filter_max_diff, prob_85 = 0.85, prob_95 = 0.95)
uncertainty3_m5v2 <- extract_summary(m5v2_e_prop_diff, prob_85 = 0.85, prob_95 = 0.95)
uncertainty4_m5v2 <- extract_summary(m5v2_R_c_diff, prob_85 = 0.85, prob_95 = 0.95)
uncertainty5_m5v2 <- extract_summary(m5v2_G_c_diff, prob_85 = 0.85, prob_95 = 0.95)
uncertainty6_m5v2 <- extract_summary(m5v2_B_c_diff, prob_85 = 0.85, prob_95 = 0.95)
Show code cell content
# Save tables
write.table(uncertainty1_m5v2, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_uncertainty_m5v2_e_max_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(uncertainty2_m5v2, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_uncertainty_m5v2_Filter_max_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(uncertainty3_m5v2, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_uncertainty_m5v2_e_prop_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(uncertainty4_m5v2, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_uncertainty_m5v2_R_c_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(uncertainty5_m5v2, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_uncertainty_m5v2_G_c_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
write.table(uncertainty6_m5v2, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_uncertainty_m5v2_B_c_diff.csv", sep = ",", row.names = TRUE, col.names = TRUE)
Show code cell source
# Convert each data frame to a plain HTML table string
table_1_uncertainty_m5v2 <- minimal_html_table(uncertainty1_m5v2, caption = "Uncertainty values - e_max Background Matching")
table_2_uncertainty_m5v2 <- minimal_html_table(uncertainty2_m5v2, caption = "Uncertainty values - Filter_max Background Matching")
table_3_uncertainty_m5v2 <- minimal_html_table(uncertainty3_m5v2, caption = "Uncertainty values - e_prop Background Matching")
table_4_uncertainty_m5v2 <- minimal_html_table(uncertainty4_m5v2, caption = "Uncertainty values - R Background Matching")
table_5_uncertainty_m5v2 <- minimal_html_table(uncertainty5_m5v2, caption = "Uncertainty values - G Background Matching")
table_6_uncertainty_m5v2 <- minimal_html_table(uncertainty6_m5v2, caption = "Uncertainty values - B Background Matching")
my_tabs_uncertainty_m5v2 <- '
<style>
/* Basic container styling */
.tabs-container {
width: 100%;
margin: 1em 0;
}
/* Hide the radio inputs (we only show their labels as tabs) */
.tabs-container input[type="radio"] {
display: none;
}
/* The “tab-label” styling: looks like a tab */
.tab-label {
display: inline-block;
padding: 10px;
margin-right: 2px;
background: #eee;
border: 1px solid #ccc;
cursor: pointer;
border-bottom: none;
}
/* The active tab label */
.tab-label-active {
background: #fff;
}
/* The panel that holds table content */
.tab-content {
border: 1px solid #ccc;
padding: 10px;
display: none;
}
/* For each radio input, show its corresponding content when checked */
#tab1_uncertainty_m5v2:checked ~ #content1_uncertainty_m5v2,
#tab2_uncertainty_m5v2:checked ~ #content2_uncertainty_m5v2,
#tab3_uncertainty_m5v2:checked ~ #content3_uncertainty_m5v2,
#tab4_uncertainty_m5v2:checked ~ #content4_uncertainty_m5v2,
#tab5_uncertainty_m5v2:checked ~ #content5_uncertainty_m5v2,
#tab6_uncertainty_m5v2:checked ~ #content6_uncertainty_m5v2 {
display: block;
}
/* Also style the label of the checked radio as “active” using the :checked + label technique */
#tab1_uncertainty_m5v2:checked + label[for="tab1_uncertainty_m5v2"],
#tab2_uncertainty_m5v2:checked + label[for="tab2_uncertainty_m5v2"],
#tab3_uncertainty_m5v2:checked + label[for="tab3_uncertainty_m5v2"],
#tab4_uncertainty_m5v2:checked + label[for="tab4_uncertainty_m5v2"],
#tab5_uncertainty_m5v2:checked + label[for="tab5_uncertainty_m5v2"],
#tab6_uncertainty_m5v2:checked + label[for="tab6_uncertainty_m5v2"] {
background: #fff;
border-bottom: none;
}
</style>
<div class="tabs-container">
<!-- 1) Tab radio + label -->
<input type="radio" name="tabs_uncertainty_m5v2" id="tab1_uncertainty_m5v2" checked>
<label class="tab-label" for="tab1_uncertainty_m5v2">Table 1</label>
<!-- 2) Tab radio + label -->
<input type="radio" name="tabs_uncertainty_m5v2" id="tab2_uncertainty_m5v2">
<label class="tab-label" for="tab2_uncertainty_m5v2">Table 2</label>
<!-- 3) Tab radio + label -->
<input type="radio" name="tabs_uncertainty_m5v2" id="tab3_uncertainty_m5v2">
<label class="tab-label" for="tab3_uncertainty_m5v2">Table 3</label>
<!-- 4) Tab radio + label -->
<input type="radio" name="tabs_uncertainty_m5v2" id="tab4_uncertainty_m5v2">
<label class="tab-label" for="tab4_uncertainty_m5v2">Table 4</label>
<!-- 5) Tab radio + label -->
<input type="radio" name="tabs_uncertainty_m5v2" id="tab5_uncertainty_m5v2">
<label class="tab-label" for="tab5_uncertainty_m5v2">Table 5</label>
<!-- 6) Tab radio + label -->
<input type="radio" name="tabs_uncertainty_m5v2" id="tab6_uncertainty_m5v2">
<label class="tab-label" for="tab6_uncertainty_m5v2">Table 6</label>
<!-- Content for each tab -->
<div class="tab-content" id="content1_uncertainty_m5v2">REPLACE_WITH_table_1_m5v2</div>
<div class="tab-content" id="content2_uncertainty_m5v2">REPLACE_WITH_table_2_m5v2</div>
<div class="tab-content" id="content3_uncertainty_m5v2">REPLACE_WITH_table_3_m5v2</div>
<div class="tab-content" id="content4_uncertainty_m5v2">REPLACE_WITH_table_4_m5v2</div>
<div class="tab-content" id="content5_uncertainty_m5v2">REPLACE_WITH_table_5_m5v2</div>
<div class="tab-content" id="content6_uncertainty_m5v2">REPLACE_WITH_table_6_m5v2</div>
</div>
'
# Now do the replacements for each table
my_tabs_uncertainty_m5v2 <- gsub("REPLACE_WITH_table_1_m5v2", table_1_uncertainty_m5v2, my_tabs_uncertainty_m5v2)
my_tabs_uncertainty_m5v2 <- gsub("REPLACE_WITH_table_2_m5v2", table_2_uncertainty_m5v2, my_tabs_uncertainty_m5v2)
my_tabs_uncertainty_m5v2 <- gsub("REPLACE_WITH_table_3_m5v2", table_3_uncertainty_m5v2, my_tabs_uncertainty_m5v2)
my_tabs_uncertainty_m5v2 <- gsub("REPLACE_WITH_table_4_m5v2", table_4_uncertainty_m5v2, my_tabs_uncertainty_m5v2)
my_tabs_uncertainty_m5v2 <- gsub("REPLACE_WITH_table_5_m5v2", table_5_uncertainty_m5v2, my_tabs_uncertainty_m5v2)
my_tabs_uncertainty_m5v2 <- gsub("REPLACE_WITH_table_6_m5v2", table_6_uncertainty_m5v2, my_tabs_uncertainty_m5v2)
IRdisplay::display_html(my_tabs_uncertainty_m5v2)
| Parameter | Estimate | Est.Error | l-85% CI | u-85% CI | l-95% CI | u-95% CI |
|---|---|---|---|---|---|---|
| Intercept | -0.775 | 0.0318 | -0.821 | -0.73 | -0.837 | -0.713 |
| SexF | 0.116 | 0.0117 | 0.099 | 0.132 | 0.0927 | 0.139 |
| Microhabitat_AssociationPresent | -0.0189 | 0.0126 | -0.037 | -0.000699 | -0.0432 | 0.00623 |
| Viewpoint10 | -0.0432 | 0.0139 | -0.0634 | -0.023 | -0.0704 | -0.0159 |
| Viewpoint20 | -0.0683 | 0.014 | -0.0884 | -0.0482 | -0.0952 | -0.0412 |
| Parameter | Estimate | Est.Error | l-85% CI | u-85% CI | l-95% CI | u-95% CI |
|---|---|---|---|---|---|---|
| Intercept | 2.63 | 0.0198 | 2.6 | 2.66 | 2.59 | 2.67 |
| SexF | 0.0353 | 0.00763 | 0.0244 | 0.0463 | 0.0205 | 0.0503 |
| Microhabitat_AssociationPresent | -0.0027 | 0.00811 | -0.0144 | 0.00908 | -0.0185 | 0.0133 |
| Viewpoint10 | -1.03 | 0.0092 | -1.04 | -1.02 | -1.05 | -1.01 |
| Viewpoint20 | -1.32 | 0.00911 | -1.33 | -1.31 | -1.34 | -1.3 |
| Parameter | Estimate | Est.Error | l-85% CI | u-85% CI | l-95% CI | u-95% CI |
|---|---|---|---|---|---|---|
| Intercept | -4.01 | 0.0206 | -4.04 | -3.98 | -4.05 | -3.97 |
| SexF | -0.252 | 0.00902 | -0.265 | -0.239 | -0.27 | -0.234 |
| Microhabitat_AssociationPresent | -0.0601 | 0.00946 | -0.0737 | -0.0465 | -0.0787 | -0.0413 |
| Viewpoint10 | 0.0822 | 0.0107 | 0.0666 | 0.0975 | 0.0611 | 0.103 |
| Viewpoint20 | 0.225 | 0.0107 | 0.209 | 0.24 | 0.203 | 0.245 |
| Parameter | Estimate | Est.Error | l-85% CI | u-85% CI | l-95% CI | u-95% CI |
|---|---|---|---|---|---|---|
| Intercept | 1.61 | 0.0375 | 1.56 | 1.67 | 1.54 | 1.68 |
| SexF | 0.182 | 0.0106 | 0.167 | 0.197 | 0.161 | 0.202 |
| Microhabitat_AssociationPresent | -0.0578 | 0.0114 | -0.0743 | -0.0413 | -0.0802 | -0.0355 |
| Viewpoint10 | 0.00596 | 0.0127 | -0.0122 | 0.0242 | -0.0189 | 0.0309 |
| Viewpoint20 | 0.00774 | 0.0125 | -0.00996 | 0.0261 | -0.0165 | 0.0325 |
| Parameter | Estimate | Est.Error | l-85% CI | u-85% CI | l-95% CI | u-95% CI |
|---|---|---|---|---|---|---|
| Intercept | 1.57 | 0.0342 | 1.52 | 1.61 | 1.5 | 1.63 |
| SexF | 0.159 | 0.0114 | 0.142 | 0.175 | 0.136 | 0.18 |
| Microhabitat_AssociationPresent | -0.0634 | 0.0123 | -0.081 | -0.0458 | -0.0875 | -0.0394 |
| Viewpoint10 | 0.00601 | 0.0136 | -0.0138 | 0.0255 | -0.0204 | 0.0326 |
| Viewpoint20 | 0.00657 | 0.0136 | -0.013 | 0.0263 | -0.0199 | 0.0335 |
| Parameter | Estimate | Est.Error | l-85% CI | u-85% CI | l-95% CI | u-95% CI |
|---|---|---|---|---|---|---|
| Intercept | 1.22 | 0.0365 | 1.17 | 1.27 | 1.15 | 1.29 |
| SexF | 0.252 | 0.0131 | 0.233 | 0.271 | 0.227 | 0.278 |
| Microhabitat_AssociationPresent | -0.0217 | 0.0143 | -0.0422 | -0.00112 | -0.0495 | 0.00628 |
| Viewpoint10 | 0.00619 | 0.0154 | -0.0161 | 0.0285 | -0.024 | 0.0364 |
| Viewpoint20 | 0.00704 | 0.0156 | -0.0155 | 0.0293 | -0.0238 | 0.0376 |
Check posterior probabilities#
Now we can evaluate our hypotheses using the posterior distributions of the model parameters. Remember, we are concerned with two things:
Is there a difference in background matching between females and males? (Sex effect)
Is there a difference in background matching between microhabitats? (Microhabitat effect)
We will calculate the posterior probability of these effects.
Show code cell content
##### **Posterior Probability of Sex Effect**
draws_m5v1_e_max <- as_draws_df(m5v1_e_max_diff)
draws_m5v1_Filter_max <- as_draws_df(m5v1_Filter_max_diff)
draws_m5v1_e_prop <- as_draws_df(m5v1_e_prop_diff)
draws_m5v1_R_c <- as_draws_df(m5v1_R_c_diff)
draws_m5v1_G_c <- as_draws_df(m5v1_G_c_diff)
draws_m5v1_B_c <- as_draws_df(m5v1_B_c_diff)
draws_m5v2_e_max <- as_draws_df(m5v2_e_max_diff)
draws_m5v2_Filter_max <- as_draws_df(m5v2_Filter_max_diff)
draws_m5v2_e_prop <- as_draws_df(m5v2_e_prop_diff)
draws_m5v2_R_c <- as_draws_df(m5v2_R_c_diff)
draws_m5v2_G_c <- as_draws_df(m5v2_G_c_diff)
draws_m5v2_B_c <- as_draws_df(m5v2_B_c_diff)
# Posterior probability that difference between pod and background is a **higher** for females than males (reference sex)
pp_m5v1_e_max_positive <- mean(draws_m5v1_e_max$b_SexF > 0)
pp_m5v1_Filter_max_positive <- mean(draws_m5v1_Filter_max$b_SexF > 0)
pp_m5v1_e_prop_positive <- mean(draws_m5v1_e_prop$b_SexF > 0)
pp_m5v1_R_c_positive <- mean(draws_m5v1_R_c$b_SexF > 0)
pp_m5v1_G_c_positive <- mean(draws_m5v1_G_c$b_SexF > 0)
pp_m5v1_B_c_positive <- mean(draws_m5v1_B_c$b_SexF > 0)
pp_m5v2_e_max_positive <- mean(draws_m5v2_e_max$b_SexF > 0)
pp_m5v2_Filter_max_positive <- mean(draws_m5v2_Filter_max$b_SexF > 0)
pp_m5v2_e_prop_positive <- mean(draws_m5v2_e_prop$b_SexF > 0)
pp_m5v2_R_c_positive <- mean(draws_m5v2_R_c$b_SexF > 0)
pp_m5v2_G_c_positive <- mean(draws_m5v2_G_c$b_SexF > 0)
pp_m5v2_B_c_positive <- mean(draws_m5v2_B_c$b_SexF > 0)
# Posterior probability that difference between pod and background is a **lower** for females than males (reference sex)
pp_m5v1_e_max_negative <- mean(draws_m5v1_e_max$b_SexF < 0)
pp_m5v1_Filter_max_negative <- mean(draws_m5v1_Filter_max$b_SexF < 0)
pp_m5v1_e_prop_negative <- mean(draws_m5v1_e_prop$b_SexF < 0)
pp_m5v1_R_c_negative <- mean(draws_m5v1_R_c$b_SexF < 0)
pp_m5v1_G_c_negative <- mean(draws_m5v1_G_c$b_SexF < 0)
pp_m5v1_B_c_negative <- mean(draws_m5v1_B_c$b_SexF < 0)
pp_m5v2_e_max_negative <- mean(draws_m5v2_e_max$b_SexF < 0)
pp_m5v2_Filter_max_negative <- mean(draws_m5v2_Filter_max$b_SexF < 0)
pp_m5v2_e_prop_negative <- mean(draws_m5v2_e_prop$b_SexF < 0)
pp_m5v2_R_c_negative <- mean(draws_m5v2_R_c$b_SexF < 0)
pp_m5v2_G_c_negative <- mean(draws_m5v2_G_c$b_SexF < 0)
pp_m5v2_B_c_negative <- mean(draws_m5v2_B_c$b_SexF < 0)
# Create a summary table
pp_summary_m5v1_sex <- data.frame(
Hypothesis = c(
"P(Sex effect e_max diff > 0)",
"P(Sex effect Filter_max diff > 0)",
"P(Sex effect e_prop diff > 0)",
"P(Sex effect R_c diff > 0)",
"P(Sex effect G_c diff > 0)",
"P(Sex effect B_c diff > 0)",
"P(Sex effect e_max diff < 0)",
"P(Sex effect Filter_max diff < 0)",
"P(Sex effect e_prop diff < 0)",
"P(Sex effect R_c diff < 0)",
"P(Sex effect G_c diff < 0)",
"P(Sex effect B_c diff < 0)"
),
PosteriorProbability = c(
pp_m5v1_e_max_positive,
pp_m5v1_Filter_max_positive,
pp_m5v1_e_prop_positive,
pp_m5v1_R_c_positive,
pp_m5v1_G_c_positive,
pp_m5v1_B_c_positive,
pp_m5v1_e_max_negative,
pp_m5v1_Filter_max_negative,
pp_m5v1_e_prop_negative,
pp_m5v1_R_c_negative,
pp_m5v1_G_c_negative,
pp_m5v1_B_c_negative
)
)
pp_summary_m5v2_sex <- data.frame(
Hypothesis = c(
"P(Sex effect e_max diff > 0)",
"P(Sex effect Filter_max diff > 0)",
"P(Sex effect e_prop diff > 0)",
"P(Sex effect R_c diff > 0)",
"P(Sex effect G_c diff > 0)",
"P(Sex effect B_c diff > 0)",
"P(Sex effect e_max diff < 0)",
"P(Sex effect Filter_max diff < 0)",
"P(Sex effect e_prop diff < 0)",
"P(Sex effect R_c diff < 0)",
"P(Sex effect G_c diff < 0)",
"P(Sex effect B_c diff < 0)"
),
PosteriorProbability = c(
pp_m5v2_e_max_positive,
pp_m5v2_Filter_max_positive,
pp_m5v2_e_prop_positive,
pp_m5v2_R_c_positive,
pp_m5v2_G_c_positive,
pp_m5v2_B_c_positive,
pp_m5v2_e_max_negative,
pp_m5v2_Filter_max_negative,
pp_m5v2_e_prop_negative,
pp_m5v2_R_c_negative,
pp_m5v2_G_c_negative,
pp_m5v2_B_c_negative
)
)
# Save combined table
pp_summary_sex_combined <- rbind(pp_summary_m5v1_sex, pp_summary_m5v2_sex)
write.table(pp_summary_sex_combined, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_postprob_m5_sex.csv", sep = ",", row.names = TRUE, col.names = TRUE)
Show code cell source
# Convert each to HTML
html_m5v1_sex <- minimal_html_table(pp_summary_m5v1_sex, caption = "Posterior probabilities (Sex)")
html_m5v2_sex <- minimal_html_table(pp_summary_m5v2_sex, caption = "Posterior probabilities (Sex)")
# Tabs layout
my_tabs_pp_sex <- '
<style>
.tabs-container {
width: 100%;
margin: 1em 0;
}
.tabs-container input[type="radio"] {
display: none;
}
.tab-label {
display: inline-block;
padding: 10px;
margin-right: 2px;
background: #eee;
border: 1px solid #ccc;
cursor: pointer;
border-bottom: none;
}
.tab-content {
border: 1px solid #ccc;
padding: 10px;
display: none;
}
#tab_m5v1_sex:checked ~ #content_m5v1_sex,
#tab_m5v2_sex:checked ~ #content_m5v2_sex {
display: block;
}
#tab_m5v1_sex:checked + label[for="tab_m5v1_sex"],
#tab_m5v2_sex:checked + label[for="tab_m5v2_sex"] {
background: #fff;
border-bottom: none;
}
</style>
<div class="tabs-container">
<input type="radio" name="tabs_pp_sex" id="tab_m5v1_sex" checked>
<label class="tab-label" for="tab_m5v1_sex">Microhabitat Model</label>
<input type="radio" name="tabs_pp_sex" id="tab_m5v2_sex">
<label class="tab-label" for="tab_m5v2_sex">Microhabitat_Association Model</label>
<div class="tab-content" id="content_m5v1_sex">REPLACE_WITH_m5v1_sex</div>
<div class="tab-content" id="content_m5v2_sex">REPLACE_WITH_m5v2_sex</div>
</div>
'
# Inject tables
my_tabs_pp_sex <- gsub("REPLACE_WITH_m5v1_sex", html_m5v1_sex, my_tabs_pp_sex)
my_tabs_pp_sex <- gsub("REPLACE_WITH_m5v2_sex", html_m5v2_sex, my_tabs_pp_sex)
# Display
IRdisplay::display_html(my_tabs_pp_sex)
| Hypothesis | PosteriorProbability |
|---|---|
| P(Sex effect e_max diff > 0) | 1 |
| P(Sex effect Filter_max diff > 0) | 1 |
| P(Sex effect e_prop diff > 0) | 0 |
| P(Sex effect R_c diff > 0) | 1 |
| P(Sex effect G_c diff > 0) | 1 |
| P(Sex effect B_c diff > 0) | 1 |
| P(Sex effect e_max diff < 0) | 0 |
| P(Sex effect Filter_max diff < 0) | 0 |
| P(Sex effect e_prop diff < 0) | 1 |
| P(Sex effect R_c diff < 0) | 0 |
| P(Sex effect G_c diff < 0) | 0 |
| P(Sex effect B_c diff < 0) | 0 |
| Hypothesis | PosteriorProbability |
|---|---|
| P(Sex effect e_max diff > 0) | 1 |
| P(Sex effect Filter_max diff > 0) | 1 |
| P(Sex effect e_prop diff > 0) | 0 |
| P(Sex effect R_c diff > 0) | 1 |
| P(Sex effect G_c diff > 0) | 1 |
| P(Sex effect B_c diff > 0) | 1 |
| P(Sex effect e_max diff < 0) | 0 |
| P(Sex effect Filter_max diff < 0) | 0 |
| P(Sex effect e_prop diff < 0) | 1 |
| P(Sex effect R_c diff < 0) | 0 |
| P(Sex effect G_c diff < 0) | 0 |
| P(Sex effect B_c diff < 0) | 0 |
Show code cell content
##### **Posterior Probability of Microhabitat Effect**
draws_m5v1_e_max <- as_draws_df(m5v1_e_max_diff)
draws_m5v1_Filter_max <- as_draws_df(m5v1_Filter_max_diff)
draws_m5v1_e_prop <- as_draws_df(m5v1_e_prop_diff)
draws_m5v1_R_c <- as_draws_df(m5v1_R_c_diff)
draws_m5v1_G_c <- as_draws_df(m5v1_G_c_diff)
draws_m5v1_B_c <- as_draws_df(m5v1_B_c_diff)
draws_m5v2_e_max <- as_draws_df(m5v2_e_max_diff)
draws_m5v2_Filter_max <- as_draws_df(m5v2_Filter_max_diff)
draws_m5v2_e_prop <- as_draws_df(m5v2_e_prop_diff)
draws_m5v2_R_c <- as_draws_df(m5v2_R_c_diff)
draws_m5v2_G_c <- as_draws_df(m5v2_G_c_diff)
draws_m5v2_B_c <- as_draws_df(m5v2_B_c_diff)
# Posterior probability that difference between pod and background is a **higher** for females than males (reference sex)
pp_RedAlgae_m5v1_e_max_positive <- mean(draws_m5v1_e_max$b_MicrohabitatRed_Algae > 0)
pp_RedAlgae_m5v1_Filter_max_positive <- mean(draws_m5v1_Filter_max$b_MicrohabitatRed_Algae > 0)
pp_RedAlgae_m5v1_e_prop_positive <- mean(draws_m5v1_e_prop$b_MicrohabitatRed_Algae > 0)
pp_RedAlgae_m5v1_R_c_positive <- mean(draws_m5v1_R_c$b_MicrohabitatRed_Algae > 0)
pp_RedAlgae_m5v1_G_c_positive <- mean(draws_m5v1_G_c$b_MicrohabitatRed_Algae > 0)
pp_RedAlgae_m5v1_B_c_positive <- mean(draws_m5v1_B_c$b_MicrohabitatRed_Algae > 0)
pp_HydroidBryozoa_m5v1_e_max_positive <- mean(draws_m5v1_e_max$b_MicrohabitatHydroid_Bryozoa > 0)
pp_HydroidBryozoa_m5v1_Filter_max_positive <- mean(draws_m5v1_Filter_max$b_MicrohabitatHydroid_Bryozoa > 0)
pp_HydroidBryozoa_m5v1_e_prop_positive <- mean(draws_m5v1_e_prop$b_MicrohabitatHydroid_Bryozoa > 0)
pp_HydroidBryozoa_m5v1_R_c_positive <- mean(draws_m5v1_R_c$b_MicrohabitatHydroid_Bryozoa > 0)
pp_HydroidBryozoa_m5v1_G_c_positive <- mean(draws_m5v1_G_c$b_MicrohabitatHydroid_Bryozoa > 0)
pp_HydroidBryozoa_m5v1_B_c_positive <- mean(draws_m5v1_B_c$b_MicrohabitatHydroid_Bryozoa > 0)
pp_Microhabitat_AssociationPresent_m5v2_e_max_positive <- mean(draws_m5v2_e_max$b_Microhabitat_AssociationPresent > 0)
pp_Microhabitat_AssociationPresent_m5v2_Filter_max_positive <- mean(draws_m5v2_Filter_max$b_Microhabitat_AssociationPresent > 0)
pp_Microhabitat_AssociationPresent_m5v2_e_prop_positive <- mean(draws_m5v2_e_prop$b_Microhabitat_AssociationPresent > 0)
pp_Microhabitat_AssociationPresent_m5v2_R_c_positive <- mean(draws_m5v2_R_c$b_Microhabitat_AssociationPresent > 0)
pp_Microhabitat_AssociationPresent_m5v2_G_c_positive <- mean(draws_m5v2_G_c$b_Microhabitat_AssociationPresent > 0)
pp_Microhabitat_AssociationPresent_m5v2_B_c_positive <- mean(draws_m5v2_B_c$b_Microhabitat_AssociationPresent > 0)
# Posterior probability that difference between pod and background is a **lower** for females than males (reference sex)
pp_RedAlgae_m5v1_e_max_negative <- mean(draws_m5v1_e_max$b_MicrohabitatRed_Algae < 0)
pp_RedAlgae_m5v1_Filter_max_negative <- mean(draws_m5v1_Filter_max$b_MicrohabitatRed_Algae < 0)
pp_RedAlgae_m5v1_e_prop_negative <- mean(draws_m5v1_e_prop$b_MicrohabitatRed_Algae < 0)
pp_RedAlgae_m5v1_R_c_negative <- mean(draws_m5v1_R_c$b_MicrohabitatRed_Algae < 0)
pp_RedAlgae_m5v1_G_c_negative <- mean(draws_m5v1_G_c$b_MicrohabitatRed_Algae < 0)
pp_RedAlgae_m5v1_B_c_negative <- mean(draws_m5v1_B_c$b_MicrohabitatRed_Algae < 0)
pp_HydroidBryozoa_m5v1_e_max_negative <- mean(draws_m5v1_e_max$b_MicrohabitatHydroid_Bryozoa < 0)
pp_HydroidBryozoa_m5v1_Filter_max_negative <- mean(draws_m5v1_Filter_max$b_MicrohabitatHydroid_Bryozoa < 0)
pp_HydroidBryozoa_m5v1_e_prop_negative <- mean(draws_m5v1_e_prop$b_MicrohabitatHydroid_Bryozoa < 0)
pp_HydroidBryozoa_m5v1_R_c_negative <- mean(draws_m5v1_R_c$b_MicrohabitatHydroid_Bryozoa < 0)
pp_HydroidBryozoa_m5v1_G_c_negative <- mean(draws_m5v1_G_c$b_MicrohabitatHydroid_Bryozoa < 0)
pp_HydroidBryozoa_m5v1_B_c_negative <- mean(draws_m5v1_B_c$b_MicrohabitatHydroid_Bryozoa < 0)
pp_Microhabitat_AssociationPresent_m5v2_e_max_negative <- mean(draws_m5v2_e_max$b_Microhabitat_AssociationPresent < 0)
pp_Microhabitat_AssociationPresent_m5v2_Filter_max_negative <- mean(draws_m5v2_Filter_max$b_Microhabitat_AssociationPresent < 0)
pp_Microhabitat_AssociationPresent_m5v2_e_prop_negative <- mean(draws_m5v2_e_prop$b_Microhabitat_AssociationPresent < 0)
pp_Microhabitat_AssociationPresent_m5v2_R_c_negative <- mean(draws_m5v2_R_c$b_Microhabitat_AssociationPresent < 0)
pp_Microhabitat_AssociationPresent_m5v2_G_c_negative <- mean(draws_m5v2_G_c$b_Microhabitat_AssociationPresent < 0)
pp_Microhabitat_AssociationPresent_m5v2_B_c_negative <- mean(draws_m5v2_B_c$b_Microhabitat_AssociationPresent < 0)
# Create a summary table
pp_summary_m5v1_microhabitat <- data.frame(
Hypothesis = c(
"P(Microhabitat Red Algae effect e_max diff > 0)",
"P(Microhabitat Red Algae effect Filter_max diff > 0)",
"P(Microhabitat Red Algae effect e_prop diff > 0)",
"P(Microhabitat Red Algae effect R_c diff > 0)",
"P(Microhabitat Red Algae effect G_c diff > 0)",
"P(Microhabitat Red Algae effect B_c diff > 0)",
"P(Microhabitat Hydroid with Bryozoan effect e_max diff > 0)",
"P(Microhabitat Hydroid with Bryozoan effect Filter_max diff > 0)",
"P(Microhabitat Hydroid with Bryozoan effect e_prop diff > 0)",
"P(Microhabitat Hydroid with Bryozoan effect R_c diff > 0)",
"P(Microhabitat Hydroid with Bryozoan effect G_c diff > 0)",
"P(Microhabitat Hydroid with Bryozoan effect B_c diff > 0)",
"P(Microhabitat Red Algae effect e_max diff < 0)",
"P(Microhabitat Red Algae effect Filter_max diff < 0)",
"P(Microhabitat Red Algae effect e_prop diff < 0)",
"P(Microhabitat Red Algae effect R_c diff < 0)",
"P(Microhabitat Red Algae effect G_c diff < 0)",
"P(Microhabitat Red Algae effect B_c diff < 0)",
"P(Microhabitat Hydroid with Bryozoan effect e_max diff < 0)",
"P(Microhabitat Hydroid with Bryozoan effect Filter_max diff < 0)",
"P(Microhabitat Hydroid with Bryozoan effect e_prop diff < 0)",
"P(Microhabitat Hydroid with Bryozoan effect R_c diff < 0)",
"P(Microhabitat Hydroid with Bryozoan effect G_c diff < 0)",
"P(Microhabitat Hydroid with Bryozoan effect B_c diff < 0)"
),
PosteriorProbability = c(
pp_RedAlgae_m5v1_e_max_positive,
pp_RedAlgae_m5v1_Filter_max_positive,
pp_RedAlgae_m5v1_e_prop_positive,
pp_RedAlgae_m5v1_R_c_positive,
pp_RedAlgae_m5v1_G_c_positive,
pp_RedAlgae_m5v1_B_c_positive,
pp_HydroidBryozoa_m5v1_e_max_positive,
pp_HydroidBryozoa_m5v1_Filter_max_positive,
pp_HydroidBryozoa_m5v1_e_prop_positive,
pp_HydroidBryozoa_m5v1_R_c_positive,
pp_HydroidBryozoa_m5v1_G_c_positive,
pp_HydroidBryozoa_m5v1_B_c_positive,
pp_RedAlgae_m5v1_e_max_negative,
pp_RedAlgae_m5v1_Filter_max_negative,
pp_RedAlgae_m5v1_e_prop_negative,
pp_RedAlgae_m5v1_R_c_negative,
pp_RedAlgae_m5v1_G_c_negative,
pp_RedAlgae_m5v1_B_c_negative,
pp_HydroidBryozoa_m5v1_e_max_negative,
pp_HydroidBryozoa_m5v1_Filter_max_negative,
pp_HydroidBryozoa_m5v1_e_prop_negative,
pp_HydroidBryozoa_m5v1_R_c_negative,
pp_HydroidBryozoa_m5v1_G_c_negative,
pp_HydroidBryozoa_m5v1_B_c_negative
)
)
pp_summary_m5v2_microhabitat <- data.frame(
Hypothesis = c(
"P(Microhabitat Association Present effect e_max diff > 0)",
"P(Microhabitat Association Present effect Filter_max diff > 0)",
"P(Microhabitat Association Present effect e_prop diff > 0)",
"P(Microhabitat Association Present effect R_c diff > 0)",
"P(Microhabitat Association Present effect G_c diff > 0)",
"P(Microhabitat Association Present effect B_c diff > 0)",
"P(Microhabitat Association Present effect e_max diff < 0)",
"P(Microhabitat Association Present effect Filter_max diff < 0)",
"P(Microhabitat Association Present effect e_prop diff < 0)",
"P(Microhabitat Association Present effect R_c diff < 0)",
"P(Microhabitat Association Present effect G_c diff < 0)",
"P(Microhabitat Association Present effect B_c diff < 0)"
),
PosteriorProbability = c(
pp_Microhabitat_AssociationPresent_m5v2_e_max_positive,
pp_Microhabitat_AssociationPresent_m5v2_Filter_max_positive,
pp_Microhabitat_AssociationPresent_m5v2_e_prop_positive,
pp_Microhabitat_AssociationPresent_m5v2_R_c_positive,
pp_Microhabitat_AssociationPresent_m5v2_G_c_positive,
pp_Microhabitat_AssociationPresent_m5v2_B_c_positive,
pp_Microhabitat_AssociationPresent_m5v2_e_max_negative,
pp_Microhabitat_AssociationPresent_m5v2_Filter_max_negative,
pp_Microhabitat_AssociationPresent_m5v2_e_prop_negative,
pp_Microhabitat_AssociationPresent_m5v2_R_c_negative,
pp_Microhabitat_AssociationPresent_m5v2_G_c_negative,
pp_Microhabitat_AssociationPresent_m5v2_B_c_negative
)
)
# Save combined table
pp_summary_sex_combined <- rbind(pp_summary_m5v1_microhabitat, pp_summary_m5v2_microhabitat)
write.table(pp_summary_sex_combined, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_postprob_m5_microhabitat.csv", sep = ",", row.names = TRUE, col.names = TRUE)
Show code cell source
# Convert each to HTML
html_m5v1_microhabitat <- minimal_html_table(pp_summary_m5v1_microhabitat, caption = "Posterior probabilities (Microhabitat)")
html_m5v2_microhabitat <- minimal_html_table(pp_summary_m5v2_microhabitat, caption = "Posterior probabilities (Microhabitat Association)")
# Tabs layout
my_tabs_pp_microhabitat <- '
<style>
.tabs-container {
width: 100%;
margin: 1em 0;
}
.tabs-container input[type="radio"] {
display: none;
}
.tab-label {
display: inline-block;
padding: 10px;
margin-right: 2px;
background: #eee;
border: 1px solid #ccc;
cursor: pointer;
border-bottom: none;
}
.tab-content {
border: 1px solid #ccc;
padding: 10px;
display: none;
}
#tab_m5v1_microhabitat:checked ~ #content_m5v1_microhabitat,
#tab_m5v2_microhabitat:checked ~ #content_m5v2_microhabitat {
display: block;
}
#tab_m5v1_microhabitat:checked + label[for="tab_m5v1_microhabitat"],
#tab_m5v2_microhabitat:checked + label[for="tab_m5v2_microhabitat"] {
background: #fff;
border-bottom: none;
}
</style>
<div class="tabs-container">
<input type="radio" name="tabs_pp_microhabitat" id="tab_m5v1_microhabitat" checked>
<label class="tab-label" for="tab_m5v1_microhabitat">Microhabitat Model</label>
<input type="radio" name="tabs_pp_microhabitat" id="tab_m5v2_microhabitat">
<label class="tab-label" for="tab_m5v2_microhabitat">Microhabitat_Association Model</label>
<div class="tab-content" id="content_m5v1_microhabitat">REPLACE_WITH_m5v1_microhabitat</div>
<div class="tab-content" id="content_m5v2_microhabitat">REPLACE_WITH_m5v2_microhabitat</div>
</div>
'
# Inject tables
my_tabs_pp_microhabitat <- gsub("REPLACE_WITH_m5v1_microhabitat", html_m5v1_microhabitat, my_tabs_pp_microhabitat)
my_tabs_pp_microhabitat <- gsub("REPLACE_WITH_m5v2_microhabitat", html_m5v2_microhabitat, my_tabs_pp_microhabitat)
# Display
IRdisplay::display_html(my_tabs_pp_microhabitat)
| Hypothesis | PosteriorProbability |
|---|---|
| P(Microhabitat Red Algae effect e_max diff > 0) | 0 |
| P(Microhabitat Red Algae effect Filter_max diff > 0) | 1 |
| P(Microhabitat Red Algae effect e_prop diff > 0) | 1 |
| P(Microhabitat Red Algae effect R_c diff > 0) | 0 |
| P(Microhabitat Red Algae effect G_c diff > 0) | 0 |
| P(Microhabitat Red Algae effect B_c diff > 0) | 0 |
| P(Microhabitat Hydroid with Bryozoan effect e_max diff > 0) | 1 |
| P(Microhabitat Hydroid with Bryozoan effect Filter_max diff > 0) | 0 |
| P(Microhabitat Hydroid with Bryozoan effect e_prop diff > 0) | 0.00035 |
| P(Microhabitat Hydroid with Bryozoan effect R_c diff > 0) | 1 |
| P(Microhabitat Hydroid with Bryozoan effect G_c diff > 0) | 1 |
| P(Microhabitat Hydroid with Bryozoan effect B_c diff > 0) | 1 |
| P(Microhabitat Red Algae effect e_max diff < 0) | 1 |
| P(Microhabitat Red Algae effect Filter_max diff < 0) | 0 |
| P(Microhabitat Red Algae effect e_prop diff < 0) | 0 |
| P(Microhabitat Red Algae effect R_c diff < 0) | 1 |
| P(Microhabitat Red Algae effect G_c diff < 0) | 1 |
| P(Microhabitat Red Algae effect B_c diff < 0) | 1 |
| P(Microhabitat Hydroid with Bryozoan effect e_max diff < 0) | 0 |
| P(Microhabitat Hydroid with Bryozoan effect Filter_max diff < 0) | 1 |
| P(Microhabitat Hydroid with Bryozoan effect e_prop diff < 0) | 0.99965 |
| P(Microhabitat Hydroid with Bryozoan effect R_c diff < 0) | 0 |
| P(Microhabitat Hydroid with Bryozoan effect G_c diff < 0) | 0 |
| P(Microhabitat Hydroid with Bryozoan effect B_c diff < 0) | 0 |
| Hypothesis | PosteriorProbability |
|---|---|
| P(Microhabitat Association Present effect e_max diff > 0) | 0.0683 |
| P(Microhabitat Association Present effect Filter_max diff > 0) | 0.3675 |
| P(Microhabitat Association Present effect e_prop diff > 0) | 0 |
| P(Microhabitat Association Present effect R_c diff > 0) | 0 |
| P(Microhabitat Association Present effect G_c diff > 0) | 0 |
| P(Microhabitat Association Present effect B_c diff > 0) | 0.0646 |
| P(Microhabitat Association Present effect e_max diff < 0) | 0.9317 |
| P(Microhabitat Association Present effect Filter_max diff < 0) | 0.6325 |
| P(Microhabitat Association Present effect e_prop diff < 0) | 1 |
| P(Microhabitat Association Present effect R_c diff < 0) | 1 |
| P(Microhabitat Association Present effect G_c diff < 0) | 1 |
| P(Microhabitat Association Present effect B_c diff < 0) | 0.9354 |
Show code cell content
##### **Posterior Probability of Viewpoint Effect**
draws_m5v1_e_max <- as_draws_df(m5v1_e_max_diff)
draws_m5v1_Filter_max <- as_draws_df(m5v1_Filter_max_diff)
draws_m5v1_e_prop <- as_draws_df(m5v1_e_prop_diff)
draws_m5v1_R_c <- as_draws_df(m5v1_R_c_diff)
draws_m5v1_G_c <- as_draws_df(m5v1_G_c_diff)
draws_m5v1_B_c <- as_draws_df(m5v1_B_c_diff)
draws_m5v2_e_max <- as_draws_df(m5v2_e_max_diff)
draws_m5v2_Filter_max <- as_draws_df(m5v2_Filter_max_diff)
draws_m5v2_e_prop <- as_draws_df(m5v2_e_prop_diff)
draws_m5v2_R_c <- as_draws_df(m5v2_R_c_diff)
draws_m5v2_G_c <- as_draws_df(m5v2_G_c_diff)
draws_m5v2_B_c <- as_draws_df(m5v2_B_c_diff)
# Posterior probability that difference between pod and background is a **higher** for females than males (reference sex)
pp_Viewpoint10_m5v1_e_max_positive <- mean(draws_m5v1_e_max$b_Viewpoint10 > 0)
pp_Viewpoint10_m5v1_Filter_max_positive <- mean(draws_m5v1_Filter_max$b_Viewpoint10 > 0)
pp_Viewpoint10_m5v1_e_prop_positive <- mean(draws_m5v1_e_prop$b_Viewpoint10 > 0)
pp_Viewpoint10_m5v1_R_c_positive <- mean(draws_m5v1_R_c$b_Viewpoint10 > 0)
pp_Viewpoint10_m5v1_G_c_positive <- mean(draws_m5v1_G_c$b_Viewpoint10 > 0)
pp_Viewpoint10_m5v1_B_c_positive <- mean(draws_m5v1_B_c$b_Viewpoint10 > 0)
pp_Viewpoint20_m5v1_e_max_positive <- mean(draws_m5v1_e_max$b_Viewpoint20 > 0)
pp_Viewpoint20_m5v1_Filter_max_positive <- mean(draws_m5v1_Filter_max$b_Viewpoint20 > 0)
pp_Viewpoint20_m5v1_e_prop_positive <- mean(draws_m5v1_e_prop$b_Viewpoint20 > 0)
pp_Viewpoint20_m5v1_R_c_positive <- mean(draws_m5v1_R_c$b_Viewpoint20 > 0)
pp_Viewpoint20_m5v1_G_c_positive <- mean(draws_m5v1_G_c$b_Viewpoint20 > 0)
pp_Viewpoint20_m5v1_B_c_positive <- mean(draws_m5v1_B_c$b_Viewpoint20 > 0)
pp_Viewpoint10_m5v2_e_max_positive <- mean(draws_m5v2_e_max$b_Viewpoint10 > 0)
pp_Viewpoint10_m5v2_Filter_max_positive <- mean(draws_m5v2_Filter_max$b_Viewpoint10 > 0)
pp_Viewpoint10_m5v2_e_prop_positive <- mean(draws_m5v2_e_prop$b_Viewpoint10 > 0)
pp_Viewpoint10_m5v2_R_c_positive <- mean(draws_m5v2_R_c$b_Viewpoint10 > 0)
pp_Viewpoint10_m5v2_G_c_positive <- mean(draws_m5v2_G_c$b_Viewpoint10 > 0)
pp_Viewpoint10_m5v2_B_c_positive <- mean(draws_m5v2_B_c$b_Viewpoint10 > 0)
pp_Viewpoint20_m5v2_e_max_positive <- mean(draws_m5v2_e_max$b_Viewpoint20 > 0)
pp_Viewpoint20_m5v2_Filter_max_positive <- mean(draws_m5v2_Filter_max$b_Viewpoint20 > 0)
pp_Viewpoint20_m5v2_e_prop_positive <- mean(draws_m5v2_e_prop$b_Viewpoint20 > 0)
pp_Viewpoint20_m5v2_R_c_positive <- mean(draws_m5v2_R_c$b_Viewpoint20 > 0)
pp_Viewpoint20_m5v2_G_c_positive <- mean(draws_m5v2_G_c$b_Viewpoint20 > 0)
pp_Viewpoint20_m5v2_B_c_positive <- mean(draws_m5v2_B_c$b_Viewpoint20 > 0)
# Posterior probability that difference between pod and background is a **lower** for females than males (reference sex)
pp_Viewpoint10_m5v1_e_max_negative <- mean(draws_m5v1_e_max$b_Viewpoint10 < 0)
pp_Viewpoint10_m5v1_Filter_max_negative <- mean(draws_m5v1_Filter_max$b_Viewpoint10 < 0)
pp_Viewpoint10_m5v1_e_prop_negative <- mean(draws_m5v1_e_prop$b_Viewpoint10 < 0)
pp_Viewpoint10_m5v1_R_c_negative <- mean(draws_m5v1_R_c$b_Viewpoint10 < 0)
pp_Viewpoint10_m5v1_G_c_negative <- mean(draws_m5v1_G_c$b_Viewpoint10 < 0)
pp_Viewpoint10_m5v1_B_c_negative <- mean(draws_m5v1_B_c$b_Viewpoint10 < 0)
pp_Viewpoint20_m5v1_e_max_negative <- mean(draws_m5v1_e_max$b_Viewpoint20 < 0)
pp_Viewpoint20_m5v1_Filter_max_negative <- mean(draws_m5v1_Filter_max$b_Viewpoint20 < 0)
pp_Viewpoint20_m5v1_e_prop_negative <- mean(draws_m5v1_e_prop$b_Viewpoint20 < 0)
pp_Viewpoint20_m5v1_R_c_negative <- mean(draws_m5v1_R_c$b_Viewpoint20 < 0)
pp_Viewpoint20_m5v1_G_c_negative <- mean(draws_m5v1_G_c$b_Viewpoint20 < 0)
pp_Viewpoint20_m5v1_B_c_negative <- mean(draws_m5v1_B_c$b_Viewpoint20 < 0)
pp_Viewpoint10_m5v2_e_max_negative <- mean(draws_m5v2_e_max$b_Viewpoint10 < 0)
pp_Viewpoint10_m5v2_Filter_max_negative <- mean(draws_m5v2_Filter_max$b_Viewpoint10 < 0)
pp_Viewpoint10_m5v2_e_prop_negative <- mean(draws_m5v2_e_prop$b_Viewpoint10 < 0)
pp_Viewpoint10_m5v2_R_c_negative <- mean(draws_m5v2_R_c$b_Viewpoint10 < 0)
pp_Viewpoint10_m5v2_G_c_negative <- mean(draws_m5v2_G_c$b_Viewpoint10 < 0)
pp_Viewpoint10_m5v2_B_c_negative <- mean(draws_m5v2_B_c$b_Viewpoint10 < 0)
pp_Viewpoint20_m5v2_e_max_negative <- mean(draws_m5v2_e_max$b_Viewpoint20 < 0)
pp_Viewpoint20_m5v2_Filter_max_negative <- mean(draws_m5v2_Filter_max$b_Viewpoint20 < 0)
pp_Viewpoint20_m5v2_e_prop_negative <- mean(draws_m5v2_e_prop$b_Viewpoint20 < 0)
pp_Viewpoint20_m5v2_R_c_negative <- mean(draws_m5v2_R_c$b_Viewpoint20 < 0)
pp_Viewpoint20_m5v2_G_c_negative <- mean(draws_m5v2_G_c$b_Viewpoint20 < 0)
pp_Viewpoint20_m5v2_B_c_negative <- mean(draws_m5v2_B_c$b_Viewpoint20 < 0)
# Create a summary table
pp_summary_m5v1_viewpoint <- data.frame(
Hypothesis = c(
"P(Viewpoint 10 effect e_max diff > 0)",
"P(Viewpoint 10 effect Filter_max diff > 0)",
"P(Viewpoint 10 effect e_prop diff > 0)",
"P(Viewpoint 10 effect R_c diff > 0)",
"P(Viewpoint 10 effect G_c diff > 0)",
"P(Viewpoint 10 effect B_c diff > 0)",
"P(Viewpoint 20 effect e_max diff > 0)",
"P(Viewpoint 20 effect Filter_max diff > 0)",
"P(Viewpoint 20 effect e_prop diff > 0)",
"P(Viewpoint 20 effect R_c diff > 0)",
"P(Viewpoint 20 effect G_c diff > 0)",
"P(Viewpoint 20 effect B_c diff > 0)",
"P(Viewpoint 10 effect e_max diff < 0)",
"P(Viewpoint 10 effect Filter_max diff < 0)",
"P(Viewpoint 10 effect e_prop diff < 0)",
"P(Viewpoint 10 effect R_c diff < 0)",
"P(Viewpoint 10 effect G_c diff < 0)",
"P(Viewpoint 10 effect B_c diff < 0)",
"P(Viewpoint 20 effect e_max diff < 0)",
"P(Viewpoint 20 effect Filter_max diff < 0)",
"P(Viewpoint 20 effect e_prop diff < 0)",
"P(Viewpoint 20 effect R_c diff < 0)",
"P(Viewpoint 20 effect G_c diff < 0)",
"P(Viewpoint 20 effect B_c diff < 0)"
),
PosteriorProbability = c(
pp_Viewpoint10_m5v1_e_max_positive,
pp_Viewpoint10_m5v1_Filter_max_positive,
pp_Viewpoint10_m5v1_e_prop_positive,
pp_Viewpoint10_m5v1_R_c_positive,
pp_Viewpoint10_m5v1_G_c_positive,
pp_Viewpoint10_m5v1_B_c_positive,
pp_Viewpoint20_m5v1_e_max_positive,
pp_Viewpoint20_m5v1_Filter_max_positive,
pp_Viewpoint20_m5v1_e_prop_positive,
pp_Viewpoint20_m5v1_R_c_positive,
pp_Viewpoint20_m5v1_G_c_positive,
pp_Viewpoint20_m5v1_B_c_positive,
pp_Viewpoint10_m5v1_e_max_negative,
pp_Viewpoint10_m5v1_Filter_max_negative,
pp_Viewpoint10_m5v1_e_prop_negative,
pp_Viewpoint10_m5v1_R_c_negative,
pp_Viewpoint10_m5v1_G_c_negative,
pp_Viewpoint10_m5v1_B_c_negative,
pp_Viewpoint20_m5v1_e_max_negative,
pp_Viewpoint20_m5v1_Filter_max_negative,
pp_Viewpoint20_m5v1_e_prop_negative,
pp_Viewpoint20_m5v1_R_c_negative,
pp_Viewpoint20_m5v1_G_c_negative,
pp_Viewpoint20_m5v1_B_c_negative
)
)
pp_summary_m5v2_viewpoint <- data.frame(
Hypothesis = c(
"P(Viewpoint 10 effect e_max diff > 0)",
"P(Viewpoint 10 effect Filter_max diff > 0)",
"P(Viewpoint 10 effect e_prop diff > 0)",
"P(Viewpoint 10 effect R_c diff > 0)",
"P(Viewpoint 10 effect G_c diff > 0)",
"P(Viewpoint 10 effect B_c diff > 0)",
"P(Viewpoint 20 effect e_max diff > 0)",
"P(Viewpoint 20 effect Filter_max diff > 0)",
"P(Viewpoint 20 effect e_prop diff > 0)",
"P(Viewpoint 20 effect R_c diff > 0)",
"P(Viewpoint 20 effect G_c diff > 0)",
"P(Viewpoint 20 effect B_c diff > 0)",
"P(Viewpoint 10 effect e_max diff < 0)",
"P(Viewpoint 10 effect Filter_max diff < 0)",
"P(Viewpoint 10 effect e_prop diff < 0)",
"P(Viewpoint 10 effect R_c diff < 0)",
"P(Viewpoint 10 effect G_c diff < 0)",
"P(Viewpoint 10 effect B_c diff < 0)",
"P(Viewpoint 20 effect e_max diff < 0)",
"P(Viewpoint 20 effect Filter_max diff < 0)",
"P(Viewpoint 20 effect e_prop diff < 0)",
"P(Viewpoint 20 effect R_c diff < 0)",
"P(Viewpoint 20 effect G_c diff < 0)",
"P(Viewpoint 20 effect B_c diff < 0)"
),
PosteriorProbability = c(
pp_Viewpoint10_m5v2_e_max_positive,
pp_Viewpoint10_m5v2_Filter_max_positive,
pp_Viewpoint10_m5v2_e_prop_positive,
pp_Viewpoint10_m5v2_R_c_positive,
pp_Viewpoint10_m5v2_G_c_positive,
pp_Viewpoint10_m5v2_B_c_positive,
pp_Viewpoint20_m5v2_e_max_positive,
pp_Viewpoint20_m5v2_Filter_max_positive,
pp_Viewpoint20_m5v2_e_prop_positive,
pp_Viewpoint20_m5v2_R_c_positive,
pp_Viewpoint20_m5v2_G_c_positive,
pp_Viewpoint20_m5v2_B_c_positive,
pp_Viewpoint10_m5v2_e_max_negative,
pp_Viewpoint10_m5v2_Filter_max_negative,
pp_Viewpoint10_m5v2_e_prop_negative,
pp_Viewpoint10_m5v2_R_c_negative,
pp_Viewpoint10_m5v2_G_c_negative,
pp_Viewpoint10_m5v2_B_c_negative,
pp_Viewpoint20_m5v2_e_max_negative,
pp_Viewpoint20_m5v2_Filter_max_negative,
pp_Viewpoint20_m5v2_e_prop_negative,
pp_Viewpoint20_m5v2_R_c_negative,
pp_Viewpoint20_m5v2_G_c_negative,
pp_Viewpoint20_m5v2_B_c_negative
)
)
# Save combined table
pp_summary_viewpoint_combined <- rbind(pp_summary_m5v1_viewpoint, pp_summary_m5v2_viewpoint)
write.table(pp_summary_viewpoint_combined, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/tables/table_postprob_m5_viewpoint.csv", sep = ",", row.names = TRUE, col.names = TRUE)
Show code cell source
# Convert each to HTML
html_m5v1_viewpoint <- minimal_html_table(pp_summary_m5v1_viewpoint, caption = "Posterior probabilities (Viewpoint)")
html_m5v2_viewpoint <- minimal_html_table(pp_summary_m5v2_viewpoint, caption = "Posterior probabilities (Viewpoint)")
# Tabs layout
my_tabs_pp_viewpoint <- '
<style>
.tabs-container {
width: 100%;
margin: 1em 0;
}
.tabs-container input[type="radio"] {
display: none;
}
.tab-label {
display: inline-block;
padding: 10px;
margin-right: 2px;
background: #eee;
border: 1px solid #ccc;
cursor: pointer;
border-bottom: none;
}
.tab-content {
border: 1px solid #ccc;
padding: 10px;
display: none;
}
#tab_m5v1_viewpoint:checked ~ #content_m5v1_viewpoint,
#tab_m5v2_viewpoint:checked ~ #content_m5v2_viewpoint {
display: block;
}
#tab_m5v1_viewpoint:checked + label[for="tab_m5v1_viewpoint"],
#tab_m5v2_viewpoint:checked + label[for="tab_m5v2_viewpoint"] {
background: #fff;
border-bottom: none;
}
</style>
<div class="tabs-container">
<input type="radio" name="tabs_pp_viewpoint" id="tab_m5v1_viewpoint" checked>
<label class="tab-label" for="tab_m5v1_viewpoint">Microhabitat Model</label>
<input type="radio" name="tabs_pp_viewpoint" id="tab_m5v2_viewpoint">
<label class="tab-label" for="tab_m5v2_viewpoint">Microhabitat_Association Model</label>
<div class="tab-content" id="content_m5v1_viewpoint">REPLACE_WITH_m5v1_viewpoint</div>
<div class="tab-content" id="content_m5v2_viewpoint">REPLACE_WITH_m5v2_viewpoint</div>
</div>
'
# Inject tables
my_tabs_pp_viewpoint <- gsub("REPLACE_WITH_m5v1_viewpoint", html_m5v1_viewpoint, my_tabs_pp_viewpoint)
my_tabs_pp_viewpoint <- gsub("REPLACE_WITH_m5v2_viewpoint", html_m5v2_viewpoint, my_tabs_pp_viewpoint)
# Display
IRdisplay::display_html(my_tabs_pp_viewpoint)
| Hypothesis | PosteriorProbability |
|---|---|
| P(Viewpoint 10 effect e_max diff > 0) | 0.00085 |
| P(Viewpoint 10 effect Filter_max diff > 0) | 0 |
| P(Viewpoint 10 effect e_prop diff > 0) | 1 |
| P(Viewpoint 10 effect R_c diff > 0) | 0.68675 |
| P(Viewpoint 10 effect G_c diff > 0) | 0.6755 |
| P(Viewpoint 10 effect B_c diff > 0) | 0.6565 |
| P(Viewpoint 20 effect e_max diff > 0) | 0 |
| P(Viewpoint 20 effect Filter_max diff > 0) | 0 |
| P(Viewpoint 20 effect e_prop diff > 0) | 1 |
| P(Viewpoint 20 effect R_c diff > 0) | 0.73455 |
| P(Viewpoint 20 effect G_c diff > 0) | 0.6921 |
| P(Viewpoint 20 effect B_c diff > 0) | 0.66935 |
| P(Viewpoint 10 effect e_max diff < 0) | 0.99915 |
| P(Viewpoint 10 effect Filter_max diff < 0) | 1 |
| P(Viewpoint 10 effect e_prop diff < 0) | 0 |
| P(Viewpoint 10 effect R_c diff < 0) | 0.31325 |
| P(Viewpoint 10 effect G_c diff < 0) | 0.3245 |
| P(Viewpoint 10 effect B_c diff < 0) | 0.3435 |
| P(Viewpoint 20 effect e_max diff < 0) | 1 |
| P(Viewpoint 20 effect Filter_max diff < 0) | 1 |
| P(Viewpoint 20 effect e_prop diff < 0) | 0 |
| P(Viewpoint 20 effect R_c diff < 0) | 0.26545 |
| P(Viewpoint 20 effect G_c diff < 0) | 0.3079 |
| P(Viewpoint 20 effect B_c diff < 0) | 0.33065 |
| Hypothesis | PosteriorProbability |
|---|---|
| P(Viewpoint 10 effect e_max diff > 0) | 0.001 |
| P(Viewpoint 10 effect Filter_max diff > 0) | 0 |
| P(Viewpoint 10 effect e_prop diff > 0) | 1 |
| P(Viewpoint 10 effect R_c diff > 0) | 0.6793 |
| P(Viewpoint 10 effect G_c diff > 0) | 0.6702 |
| P(Viewpoint 10 effect B_c diff > 0) | 0.656 |
| P(Viewpoint 20 effect e_max diff > 0) | 0 |
| P(Viewpoint 20 effect Filter_max diff > 0) | 0 |
| P(Viewpoint 20 effect e_prop diff > 0) | 1 |
| P(Viewpoint 20 effect R_c diff > 0) | 0.72925 |
| P(Viewpoint 20 effect G_c diff > 0) | 0.68515 |
| P(Viewpoint 20 effect B_c diff > 0) | 0.6746 |
| P(Viewpoint 10 effect e_max diff < 0) | 0.999 |
| P(Viewpoint 10 effect Filter_max diff < 0) | 1 |
| P(Viewpoint 10 effect e_prop diff < 0) | 0 |
| P(Viewpoint 10 effect R_c diff < 0) | 0.3207 |
| P(Viewpoint 10 effect G_c diff < 0) | 0.3298 |
| P(Viewpoint 10 effect B_c diff < 0) | 0.344 |
| P(Viewpoint 20 effect e_max diff < 0) | 1 |
| P(Viewpoint 20 effect Filter_max diff < 0) | 1 |
| P(Viewpoint 20 effect e_prop diff < 0) | 0 |
| P(Viewpoint 20 effect R_c diff < 0) | 0.27075 |
| P(Viewpoint 20 effect G_c diff < 0) | 0.31485 |
| P(Viewpoint 20 effect B_c diff < 0) | 0.3254 |
8. Summarize results#
The predictor Sex mostly does not have a significant effect on background matching (i.e., the difference between dorsal body and microhabitat color metrics) at 85% and 95% confidence intervals. Sex MIGHT have a significant effect on pattern diversity (i.e. e_prop) and R reflectance matching between dorsal bodies and microhabitats, with females showing lower pattern diversity and R reflectance matching (i.e., higheR_c differences between dorsal body and microhabitat pattern diversity and R reflectance) at an 85% confidence interval.
Make final figures#
Generate tidy figures for the Results section of our report/paper.
Show code cell content
##### FINAL PLOT FOR WNAN PAPER
# Updated predictors list with source model per group
predictors_list <- list(
list(
name = "Sex",
model = "m1",
baseline = tibble(parameter = "Male", mean = 0),
order = c("Male", "b_SexF"),
labels = c("Male" = "Male", "b_SexF" = "Female"),
regex = "b_SexF"
),
list(
name = "Microhabitat",
model = "m1",
baseline = tibble(parameter = "Hydroids", mean = 0),
order = c("Hydroids", "b_MicrohabitatRed_Algae", "b_MicrohabitatHydroid_Bryozoa"),
labels = c(
"Hydroids" = "Hydroids",
"b_MicrohabitatRed_Algae" = "Red Algae",
"b_MicrohabitatHydroid_Bryozoa" = "Bryozoa"
),
regex = "b_Microhabitat"
),
list(
name = "Microhabitat_Association",
model = "m2", # ← uses different model
baseline = tibble(parameter = "Microhabitat_Association Not Present", mean = 0),
order = c("Microhabitat_Association Not Present", "b_Microhabitat_AssociationPresent"),
labels = c(
"Microhabitat_Association Not Present" = "Unassociated",
"b_Microhabitat_AssociationPresent" = "Associated"
),
regex = "b_Microhabitat_AssociationPresent"
),
list(
name = "Viewpoint",
model = "m1",
baseline = tibble(parameter = "No correction", mean = 0),
order = c("No correction", "b_Viewpoint10", "b_Viewpoint20"),
labels = c(
"No correction" = "No correction",
"b_Viewpoint10" = "10 mm acuity",
"b_Viewpoint20" = "20 mm acuity"
),
regex = "b_Viewpoint"
)
)
# Posterior samples by model and variable
posterior_samples_all <- list(
m1 = list(
`e[max]` = posterior_samples_m5v1_e_max_diff,
`Filter[max]` = posterior_samples_m5v1_Filter_max_diff,
`e[prop]` = posterior_samples_m5v1_e_prop_diff,
`R[c]` = posterior_samples_m5v1_R_c_diff,
`G[c]` = posterior_samples_m5v1_G_c_diff,
`B[c]` = posterior_samples_m5v1_B_c_diff
),
m2 = list(
`e[max]` = posterior_samples_m5v2_e_max_diff,
`Filter[max]` = posterior_samples_m5v2_Filter_max_diff,
`e[prop]` = posterior_samples_m5v2_e_prop_diff,
`R[c]` = posterior_samples_m5v2_R_c_diff,
`G[c]` = posterior_samples_m5v2_G_c_diff,
`B[c]` = posterior_samples_m5v2_B_c_diff
)
)
# Helper: extract posterior + baseline
extract_effects <- function(posterior_df, predictor_cfg, predictor_name) {
matched_cols <- posterior_df %>%
dplyr::select(matches(predictor_cfg$regex)) %>%
colnames()
if (length(matched_cols) == 0) {
stop(paste("No columns matched regex:", predictor_cfg$regex, "for predictor:", predictor_cfg$name))
}
draws <- posterior_df %>%
dplyr::select(all_of(matched_cols)) %>%
pivot_longer(cols = everything(), names_to = "parameter", values_to = "value") %>%
mutate(
label = predictor_cfg$labels[parameter],
predictor = predictor_name,
group = predictor_cfg$name
)
baseline <- predictor_cfg$baseline %>%
mutate(
label = predictor_cfg$labels[parameter],
predictor = predictor_name,
group = predictor_cfg$name,
value = mean
)
bind_rows(draws, baseline)
}
# Build combined plot data
plot_data <- map_dfr(
names(posterior_samples_all$m1),
function(pred_name) {
map_dfr(
predictors_list,
function(predictor_cfg) {
model_id <- predictor_cfg$model
extract_effects(
posterior_df = posterior_samples_all[[model_id]][[pred_name]],
predictor_cfg = predictor_cfg,
predictor_name = pred_name
)
}
)
}
)
plot_data$label <- factor(plot_data$label, levels = c(
"Male", "Female",
"Hydroids", "Red Algae", "Bryozoa",
"Unassociated", "Associated",
"No correction", "10 mm acuity", "20 mm acuity"
))
plot_data$predictor <- factor(
plot_data$predictor,
levels = c("e[max]", "Filter[max]", "e[prop]", "R[c]", "G[c]", "B[c]"),
labels = c(
expression(Delta*e[max]),
expression(Delta*Filter[max]),
expression(Delta*e[prop]),
expression(Delta*R[c]),
expression(Delta*G[c]),
expression(Delta*B[c])
)
)
# Plot
final_plot_posteriors_m5 <- ggplot(plot_data, aes(x = value, y = label, fill = group, color = group)) +
stat_halfeye(
.width = c(0.85, 0.95),
slab_alpha = 0.4,
interval_size_range = c(1, 0), # Inner (85%) shaded, outer (95%) invisible
normalize = "groups",
slab_linewidth = 0.6
) +
geom_point(
data = filter(plot_data, !is.na(mean)),
aes(x = mean), # ✅ Match color to group
inherit.aes = TRUE,
size = 0.75,
shape = 21
) +
geom_vline(xintercept = 0, size = 0.6, linetype = 3, color = "red") +
facet_wrap(~ predictor, ncol = 6, scales = "free_x", labeller = label_parsed) +
scale_fill_manual(
values = c(
"Sex" = "purple",
"Microhabitat" = "forestgreen",
"Microhabitat_Association" = "#ff9500",
"Viewpoint" = "steelblue"
)
) +
scale_color_manual( # ✅ Add matching colors for the points
values = c(
"Sex" = "purple",
"Microhabitat" = "forestgreen",
"Microhabitat_Association" = "#ff9500",
"Viewpoint" = "steelblue"
)
) +
labs(
x = "Effect size (log scale)",
y = NULL
) +
theme_bw(base_size = 8) +
theme(
strip.text = element_text(face = "bold", size = 10),
axis.text.y = element_text(size = 8),
panel.spacing = unit(0.5, "lines"),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
legend.position = "none"
)
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/final_plot_posteriors_m5.png", plot = final_plot_posteriors_m5, width = 9, height = 4.5, units = "in", dpi = 300)
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Warning message:
"Dropping 'draws_df' class as required metadata was removed."
Show code cell source
# Convert images to base64
final_plot_posteriors_m5 <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/final_plot_posteriors_m5.png")
# Create the HTML
html_m5_posteriors_final <- paste0("
<style>
.image-row {
display: flex;
gap: 20px;
justify-content: center;
align-items: flex-start;
}
.image-row img {
max-width: 100%;
height: auto;
border: 1px solid #ccc;
}
</style>
<div class='image-row'>
<img src='", final_plot_posteriors_m5, "' alt='Final Microhabitat Posterior Plot'>
</div>
")
# Display the HTML
IRdisplay::display_html(html_m5_posteriors_final)
